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ABSTRACT 


The objective of this study is to develop micromechanical models for predicting the 
stiffness and strength properties of textile composite materials. Micromechanical analysis 
of textile composites is possible due to the presence of a repeating unit-cell or representative 
volume element. The unit-cell is assumed to span the material continuously in all three 
dimensions. On the microscale — comparable to the scale of the unit-cell — the composite is 
heterogeneous due to the presence of the reinforcing yarn and the matrix. However, on the 
macroscale — comparable to the structural scale — the composite is assumed to be 
homogeneous and orthotropic. The homogeneous composite properties are then predicted 
from the constituent material properties and the yam geometry. 

The highlight of this study is a systematic analysis of issues involved in the finite 
element based micromechanics of textile composites. The unit-cell is discretized with 
three-dimensional finite elements, and periodic boundary conditions are imposed between 
opposite end-faces of the unit-cell. Six linearly independent deformations are applied to the 
unit-cell. From the forces acting on the unit-cell for each of the six deformations, the 
composite stiffness matrix is obtained. A similar procedure is followed to determine the 
composite coefficients of thermal expansion. The numerical procedure was tested by 
applying it for simple examples, for which the results are known. The numerical results were 
also compared with existing models for textile composites. In both cases, the results compare 
very favorably. The finite element procedure is extended to compute the thermal residual 
microstresses and to estimate the initial failure envelope for the textile composite. 

An independent finite element micromechanical analysis, analogous to the above, is 
presented for thin textile composite structures with few unit-cells in the thickness direction. 
In that case, the composite is modeled as a homogeneous plate to predict the plate stiffness 
coefficients and plate coefficients of thermal expansion. It was shown the plate properties 
could not be predicted from the corresponding three-dimensional properties. 

In addition, an approximate analytical procedure is presented to estimate the 
composite thermo-elastic constants. The procedure called the Selective Averaging Method 
is based on a judicious combination of stiffness and compliance averaging. The method is 
fast and easy to implement, and suitable for parametric studies. 



CHAPTER 1 
INTRODUCTION 

1.1 Background 

The increasing demand for lightweight yet strong and stiff structures has lead to the 
development of advanced fiber reinforced composites. These materials are not only used in 
the aerospace industry but also in a variety of commercial applications like automobile, 
marine and biomedical applications. Traditionally fibrous composites are manufactured by 
laminating several layers of unidirectional fiber tapes (Fig. 1) pre-impregnated with matrix 
material. The effective properties of the composite can be controlled by changing several 
parameters like the fiber orientation in a layer, stacking sequence, fiber and matrix material 
properties and fiber volume fraction. However, the manufacture of fibrous laminated 
composites from prepregs is labor intensive. Besides, laminated composites lack 
through-the-thickness reinforcement, and hence have poor interlaminar strength and 
fracture toughness. 
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Figure 2. Examples of woven, braided and knitted textile preforms. (Source: 
Chou and Ko, 1989) 


Recent developments in textile manufacturing processes show some promise in 
overcoming the above limitations. Textile processes such as weaving, braiding and knitting 
can turn large volumes of yam into dry preforms at a faster rate, thus reducing costs and cycle 
time. The dry preforms (Fig. 2) are impregnated with appropriate matrix material and cured 
in a mold using processes such as Resin Transfer Molding (RTM). Two-dimensional woven 
and braided mats offer increased through the thickness properties due to yam interlacing. 
The mats may be stitched using Kevlar or glass threads to provide additional reinforcement 
along the thickness direction (Sharma and Sankar, 1995). Three-dimensional woven and 
braided composites provide multidirectional reinforcement, thus directly enhancing the 
strength and stiffness in the thickness direction. Unlike laminated structures 
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three-dimensional composites do not possess the weak plane of delamination, thus giving 
increased impact resistance and fracture toughness. Textile manufacturing processes in 
conjunction with resin transfer molding are also suitable for the production of intricate 
structural forms at a reduced cycle time. This allows complex shaped structures to be 
fabricated as an integral unit, thus eliminating the use of joints and fasteners. 

With the advancements in aforementioned technologies there is a need to develop 
scientific methods of predicting the performance of the composites made using the above 
processes. There are numerous variables involved in textile processes besides the choice of 
the fiber and matrix materials. This, for example, includes (a) the number of filaments in the 
yam specified by the yam linear density and (b) the yarn architecture (description of the yam 
geometry) determined by the type of weaving or braiding processes. Thus, there is a need 
for analytical/numerical models to study the effect of these variables on the textile composite 
behavior. 

Ideally a structural engineer would like to model textile composites as a 
homogeneous anisotropic material — preferably orthotropic — so that the structural 
computations can be simplified, and also the existing computer codes can be used in the 
design. This would require the prediction of the effective (macroscopic) properties of the 
composites from the constituent material (microscopic) characteristics such as yam and 
matrix properties, yam-matrix interface characteristics and the yam architecture (Fig. 3). 
This is possible if we assume that there is a representative volume element (RVE) or an 
unit-cell that repeats its self throughout the volume of the composite, which is true in the case 
of textile composites. The unit-cell can be considered as the smallest possible building block 
for the textile composite, such that the composite can be created by assembling the unit-cell 
in all three dimensions (Fig. 4). The prediction of the effective macroscopic properties from 
the constituent material characteristics is one of the aspects of the science known as 
micromechanics. The effective properties include thermo-mechanical properties like 
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stiffness, strength and coefficients of thermal expansion as well as thermal conductivities, 
electromagnetic and other transport properties. 
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Most of the early work to determine the properties of heterogeneous materials was 
restricted to particulate (spherical) inclusions which were assumed to be isotropic. For 
example, Hashin (1962) derived expressions for the bounds for the elastic moduli of 
heterogeneous materials using variational theory. A review of analytical methods for 
predicting the effective properties of particulate composites was presented by Christensen 
(1990). However, the complex geometry of the textile preforms makes such precise 
analytical modeling difficult. 

The current literature dealing with micromechanical analyses for textile composites 
can be broadly classified into three categories: mechanics of materials type models, energy 
based approach, and finite element analysis of the unit-cell. All of the above models 
recognize that there is a unit-cell in the composite material, and they attempt to model the 
material as a homogeneous orthotropic material. In the mechanics of materials type models, 
the yams are approximated as simple structural elements, e.g., beams, plates, laminates etc., 
and their deformation behavior is assumed to be governed by the corresponding structural 
constitutive relations. The kinematics is also simplified to a great extent, and a relation 
between the overall deformation of the unit-cell and the average forces are derived. The 
energy approach is similar to the previous one, except that the strain energy in the unit-cell 
is evaluated based on some assumed displacement field, which is usually an 
oversimplification of exact displacements. The elastic constants are derived by equating the 
strain energy in the approximate model to that in an idealized homogeneous composite. 
Most energy based approaches provide bounds for the homogeneous properties, and can be 
used as a check for experimental observations or other theoretical models. The third method 
is the rigorous micromechanical analysis of the unit-cell, which often requires the use of 
numerical methods such as the finite element method, and also uses the exact 
three-dimensional constitutive relations for the yam and the matrix material. 

Ishikawa and Chou (1982a, 1982b, 1983a, 1983b, 1983c) proposed three analytical 
models for thermo-elastic properties of woven fabric composites — the mosaic model, the 
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fiber undulation model and bridging model. The mosaic model treats the composite as an 
assembly of crossply laminates, and then uses lamination theory to predict the composite 
properties. The fiber undulation model in addition takes into account the effect of yarn 
continuity in the loading direction. The bridging model was specifically developed to 
estimate the elastic behavior of satin weaves. This model simulates the effect of load 
distribution and load transfer between the yams. The above three models were developed 
for two-dimensional fabric composites. These models were then extended by Yang et al. 
(1986) in the fiber inclination model to predict the elastic properties of three-dimensional 
textile composites. The fiber inclination model was used to determine the elastic properties 
of three-dimensional angle-interlock composites and braided composites by Whitney and 
Chou (1989) and Ma et al. (1986) respectively. The results were compared with test data, 
and parametric analyses were performed to study the effect of changing geometric 
parameters. Crane and Camponeschi (1986) presented an analytical model based on classical 
lamination theory to predict the extensional stiffnesses for multi-directional braided 
composites. Naik (1994) used a stiffness averaging technique to predict woven and braided 
composite properties (explained in Chapter 3). 

Numerical modeling of the unit-cell is popular due to its ability to capture the effects 
of complicated yam architectures. For instance, Yoshino and Ohtsuka (1982) performed a 
two-dimensional finite element analysis using plane strain elements to predict the stress 
distribution within a plain-weave fabric. Dasgupta et al. (1990) and Whitcomb (1991) 
analyzed the unit-cell of a plain-weave composite using three-dimensional finite elements 
to determine the effect of the yarn geometry and yarn volume fraction on the composite 
thermo-elastic constants. Foye (1993) used inhomogeneous elements called replacement 
elements to model the unit-cell. His model can be used to predict both composite stiffness 
and strength properties. Cox et al. ( 1 994) presented a three-dimensional finite element model 
using two types of elements. The yams are modeled as two-node line elements and the rest 
of the medium as eight-node solid elements. The model was used to predict failure 
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mechanisms in angle and orthogonal interlock woven composites. Raju et al. (1990) 
compiled a review of available analytical and numerical methods for modeling textile 
composites. 


1.2 Scope of Study 

The analytical models explained above are approximate — they make use of 
assumptions similar to that in classical plate theory, and the yam geometry is greatly 
simplified. Further, these models do not consider periodic boundary conditions (see Section 
2.1.1) for continuity of displacements and tractions on opposite faces of the unit-cell. Most 
of the existing numerical models are for specific yam geometries such as simple weave 
architectures. Therefore there is a need for more general models to predict the mechanical 
properties for a textile composite, particularly for composites with complex woven and 
braided geometries. The present study illustrates two such methods. The first method 
involves rigorous finite element analysis of the unit-cell, which imposes exact displacement 
and traction boundary conditions on the end-faces of the unit-cell. The second is an 
approximate analysis in which the unit-cell is discretized into a number of elements. The 
effective (macroscale) stiffness for the composite is computed by a combination of stiffness 
and compliance averaging of the element stiffness coefficients. The advantage of the 
proposed methods is their generality — they are applicable to any textile geometry provided 
that there exists a repeatable pattern (which defines the unit-cell). Both the methods can 
compute the effective elastic constants and effective coefficients of thermal expansion for 
the homogenous composite. 

The analysis for the elastic constants of the textile composite continuum assumes that 
the unit-cell is repeatable in all three dimensions. This assumption will not be valid for thin 
textile composite structures, where there will be a only a few unit-cells in the thickness 
direction. For such applications, an alternate micromechanical analysis is proposed. The 
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unit-cell is modeled as a plate (Fig. 5) and the corresponding plate stiffness coefficients and 
plate coefficients of thermal expansion are predicted. 



Chapter 2 explains the finite element procedure to compute the macroscale 
thermo-elastic constants. The procedure is detailed for a thick textile composite in which 
there are sufficient number of unit-cells in the thickness direction. The procedure is then 
extended to a thin composite. The thin composite is first modeled as a beam and a 
two-dimensional analysis is presented to compute the beam thermo-elastic coefficients. A 
similar three-dimensional analysis is also illustrated in which the thin composite is modeled 
as a plate. An analytical procedure called the Selective Averaging Method (SAM) is 
presented in Chapter 3 to estimate the macroscale thermo-elastic constants. This procedure 
utilizes a judicious combination of stiffness and compliance averaging to compute the 
composite properties. Again for a thin textile composite the plate stiffness coefficients and 
plate coefficients of thermal expansion may be computed. 
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In Chapter 4 the microstresses computed using the finite element analysis, are used 
to construct an initial failure envelope for a textile composite modeled as homogeneous 
continuum. An independent procedure using plate theory is presented to predict the failure 
envelope for a thin textile composite. Since the composite is generally fabricated at a 
temperature greater than room temperature, thermal residual microstresses are developed 
due to the mismatch in the coefficients of thermal expansion for the constituent materials. 
A method to predict these thermal microstresses is illustrated. Chapter 5 deals with the issues 
pertaining to finite element modeling of the unit-cell. Due to the complex yam architectures, 
traditional finite element modeling becomes exceedingly difficult. Various alternative 
means of finite element modeling the unit-cell are discussed with suitable two-dimensional 
examples. Finally Chapter 6 summarizes the work done in this study. The limitations of the 
work are addressed and suggestions are given for future work in this area. 



CHAPTER 2 

FINITE ELEMENT MODELS FOR THERMO-ELASTIC CONSTANTS 

In this chapter, we demonstrate finite element based micromechanical models to 
predict the effective stiffness properties and effective coefficients of thermal expansion 
(CTE’s) for a textile composite. The macroscale properties of the composite are determined 
at a scale much larger than the dimensions of the unit-cell, but comparable to the dimensions 
of the structural component. The average stresses at a point at the structural scale will be 
called the macroscale stresses or macrostresses. The actual stresses at a point at the 
continuum level will be called the microscale stresses or microstresses. To distinguish the 
macroscale deformations and stresses from their microscale counterparts — a superscript 
“M” will be used to denote the macroscale deformations and stresses. 

2.1 Unit-Cell Analysis for Three-Dimensional Elastic Constants 

In this section a procedure to determine the three-dimensional elastic constants for 
a textile composite material is described. Consider a rectangular hexahedron as the unit-cell 
of the three-dimensional textile composite. The edges of the unit-cell are assumed to be 
parallel to the coordinate axes xj, X 2 and xj, with unit-cells repeating in all three directions. 
The length of the unit-cell in thexj direction is defined as I*. On the macroscale the composite 
is assumed to be homogeneous and orthotropic and the composite behavior is characterized 
by the following constitutive relation : 
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where {o M } and {e M } are the macrostresses and macrostrains, respectively; { a c } and [C] 
are the macroscale CTE’s and orthotropic elasticity matrix to be determined; AT M is a 
uniform temperature difference throughout the unit-cell. The macroscale elasticity matrix 
is of the form: 
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The unit-cell analysis assumes that the material is homogeneous and orthotropic on the 
macroscale. The material, therefore, is subject to a uniform state of strain in the macroscopic 
sense. The macrostresses required to create such a state of strain are computed from the finite 
element model of the unit-cell. In the microscale, all unit-cells have identical stress and strain 
fields. Continuity of microstresses across the unit-cell then requires that tractions be equal 
and opposite at corresponding points on opposite faces of the unit-cell. Since the 
displacement gradients are constant for a homogeneous deformation, the displacements at 
corresponding points on opposite faces of the unit-cell differ only by a constant. 
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2.1.1 Periodic Boundary Conditions 

The periodic boundary conditions (BC’s) consist of the periodic displacement 
boundary conditions which ensure the compatibility of displacements on opposite faces of 
the unit-cell, and the periodic traction boundary conditions to enforce the continuity of 
stresses. A macroscopically homogeneous deformation can be represented as 

uf = H ijXj i,j =1,2,3 (3) 

where Hy are the displacement gradients. Then the periodic displacement boundary 
conditions to be imposed on the faces Xi=0 and Xi=L\ are 

ufL 1*2*3) ~ u t { 0,x 2 x 3 ) — H n L j 

ufx x L 2 x 3 ) - u i (x l 0,* 3 ) = H i2 L 2 (4) 

u i (x l x 2 L 3 ) — u i (x l x 2 0) — H i3 L 3 

The traction boundary conditions to be imposed on the faces X{=0 and x;==Li are 

F i (L l x 2 ,x 3 ) = - Ffi 0,x 2 ,x 3 ) 

F i(* \,L 2 ,x 3 ) = - F i (x l 0,x 3 ) (5) 

F i {x l x 2 L 3 )= -F i {x l x 2 , 0) 

The periodic displacement and traction BC’s (for a two-dimensional problem) to be imposed 
for the deformation given by ejj =HLi and £22-yi2~0 are shown in Fig. 6. 

The above periodic BC’s are imposed in the finite element model either by using 
multi-point constraint elements or by incorporating transformation equations to eliminate 
the constrained displacements (Cook et al., 1989). These two methods to impose the periodic 
BC’s are discussed in Section 5.1.3. Both the methods require a finite element model with 
corresponding nodes on opposite faces of the unit-cell. 
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2.1.2 Determination of Three-Dimensional Elastic Constants and CTE’s 

The unit-cell is discretized with three-dimensional finite elements such that opposite 
faces of the unit-cell have identical nodes. Periodic displacement and traction boundary 
conditions are enforced between the corresponding nodes. The periodic displacement BC’s 
are imposed such that only one of the macrostrain components is non-zero; and the uniform 
temperature difference A T M is set to zero. Then, the difference in displacements between 
corresponding points on opposite faces of the unit-cell will be equal to that in a homogenous 
continuum subject to the same deformation. The average stresses (macrostresses) required 
to create such a deformation are obtained from the finite element results. Substituting the 
macrostresses and macrostrains in the composite constitutive relation (Eqn. 1) the stiffness 
coefficients in the column corresponding to the non-zero strain can be evaluated. This 
procedure is repeated for other macrostrain components (keeping the temperature difference 
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zero) to obtain the entire stiffness matrix [C]. The orthotropic elastic constants of the 
composite material can be easily determined by inverting the stiffness matrix, and 
comparing the compliance coefficients with that of an orthotropic material. 

To compute the six CTE’s, a finite temperature change T 0 is applied to all the 
elements in the unit-cell; and periodic displacement BC’s are imposed such that all the 
macrostrain components are zero. Then the composite constitutive relation Eqn. (1) will 
reduce to 

\a M ) = - [C](a c J T 0 (6) 

The macrostresses for such a deformation are computed as described below. Knowing the 
stiffness matrix [C], the composite CTE’s are found as 

M = — -f- [C] -1 (<7 m ) (7) 

1 o 


Table 1 presents the non-zero displacement BC’s imposed on the unit-cell to obtain [C] and 
the CTE’s {a c }. 

Table 1. Non-zero displacement BC’s to obtain 3-D elastic constants and CTE’s. 


stiffness coefficients to be obtained 

non-zero displacement BC’s 

first column of [C] (£ji m = 1 ) 

uj(L], x 2 , x 3 ) - uj( 0 , x 2 , x 3 ) = Lj 

second column of [C] (£ 22 M = 1 ) 

u 2 (xi, L 2 , x 3 ) - u 2 (x 1, 0 , x 3 ) = L 2 

third column of [C] (£33 M = 1 ) 

U 3(x], x 2 , L3) - U3(xj, x 2 , 0 ) = Lj 

fourth column of [C] (y 2 3 M = i) 

U 2 (xj, x 2 , L3) - U 2 (X], x 2 , 0 ) = L3/2 
U3(x], L 2 , x 3 ) - U3(xj, 0 , X3) = L 2 f 2 

fifth column of [C] (yi3 M = 1 ) 

U](X], x 2 , L 3 ) - uj(x], x 2 , 0 ) = L3/2 
u 3 (L], x 2 , x 3 ) - 1*3(0, x 2 , x 3 ) = L ]/2 

sixth column of [C] {yj 2 M = 1 ) 

U](X], L 2 , x 3 ) - Uj(xj, 0 , x 3 ) ~ L 2 /2 
u 2 (L], x 2 , X3) - u 2 ( 0 , x 2 , X3) = Lj /2 

CTE’s (A T m = 1 ) 

AT = 1 
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The macrostresses for a given deformation state can be found by the following two 
methods. In the first method, the macrostresses are obtained by averaging the nodal forces 

on each face of the unit cell. For example, the macroscale <7^ can be obtained as 

°n ~ L L y F^(L x , x 2, x t) (8) 

n 

where is the nodal force in the xj direction at the nih. node, and ^ denotes summation 

n 

over all nodes on the face xj = L j . Alternatively, the macrostresses can be computed by 
volume- averaging the corresponding microstress component in the unit-cell. Then the 
macrostress component cr^ is obtained as 

°U = yj o u (x,y,z)dV (9) 

v 

where V is the volume of the unit-cell. The microstresses are computed at the quadrature 
points, and numerically integrated over the volume in each element of the unit-cell. 


2J.3 Results and Discussion for 3-D Elastic Constants and CTE’s 


The above procedure was demonstrated for the following material systems: 
Example 1. isotropic material 

Example 2. bimaterial medium — both materials are assumed isotropic (Fig. 7) 
Example 3 . unidirectional composite with identical Poisson ratios for fiber and matrix — 

fiber and matrix materials are isotropic (Fig. 8) 

Example 4. unidirectional composite with different Poisson ratios for fiber and matrix — 

fiber and matrix materials are isotropic 

Example 5. plain-weave textile composite (Fig. 9) — yam geometry and properties 
obtained from Dasgupta et al. (1990) 

Example 6. plain-weave textile composite — yam geometry and properties obtained from 


Naik (1994) 
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Example 7. 5-hamess satin weave (Fig. 10) — yam geometry and properties obtained 
from Naik (1994) 

For the textile composite examples (examples 5-7), the yarn is assumed to be transversely 
isotropic and the matrix material is assumed isotropic. The constituent material properties 
for the examples are listed in Table 2. 


Table 2. Properties of constituent materials for examples 1-7. 


Example 1 

unit-cell size: 

E = 10 GPa, v-0.3, a=10xl0~ 6 /°C | 

0. 500x0.500x0. 256 mm 

Example 2 

layer 1 (E-glass): 
layer 2 (epoxy): 
unit-cell size: 

E]=70 GPa, v } =0.200, aj= 5x10^ /°C, V 1 =0.5 ! 

E 2 =3.50 GPa, v 2 =0.350, a 2 = 60x1 0* 6 /°C, V 2 =0.5 
0. 500x0.500x0. 256 mm 

Example 3 

fiber: 
matrix: 
unit-cell size: 

Ef=100 GPa, Vf=0.300, a/= 10x10^ /°C, Vf=0.6 
E m =10 GPa, v m =0.300, a m = lOOxlO- 6 /°C 
10x10x10 ^im 

Example 4 

fiber (E-glass): 
matrix (epoxy): 
unit-cell size: 

Ef=70 GPa, Vf=0.200, aj= 5x10^ /°C, Vf=0.6 
E m =3.50 GPa, v m =0.350, a m = 60x10^ /°C 
10x10x10 fim 

Example 5 

yarn (glass-epoxy): 

E L =58.61 GPa, E T =14.49 GPa, G LT =5.38 GPa, v LT =0.250 
VTr=0.247, a^.HxlO-t/oC, a T =22.64xlO- 6 /°C, Vf=0.26 

matrix (epoxy): 

E=3.45 GPa, v=0.37, a=69xlO~ 6 /°C 
unit-cell size: 1.680x1.680x0.228 mm 

Examples 6, 7 

yarn (graphite-epoxy): 

E L = 144.80 GPa, E T = 11.73 GPa, G LT =5.52 GPa, v LT =0.230 
vjr=0.300, a L =-0.324xl0~ 6 /°C, a T =l 4.00x1 O^C, Vf=0.64 
matrix (epoxy): 

E=3.45 GPa, v=0.35, a=40xl0~ 6 /°C 

unit-cell size: 2.822x2.822x0.2557 mm (Example 6) 

7.055x7.055x0.2557 mm (Example 7) 


Vi stands for the volume fraction of the constituent material 
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A 3-D finite element code called juTEx - 10 (pronounced as microtech) was written 
and implemented for the seven examples to compute the homogeneous elastic constants and 
CTE’s. The user manual for^TEft-lO is provided in Appendix-A. The unit-cell was assumed 
to be a rectangular hexahedron with edges along the x y- and z- axes. The unit-cell was 
divided into uniform-sized eight-node brick elements as shown in Figure 11. The number 
of divisions along the x-, y- and z- axes are specified by the user. The elements were in 
general, inhomogeneous elements, i.e., consisted of more than one constituent material. 
While computing the element stiffness matrix, the material property at the Gauss integration 
points was determined, and the corresponding elasticity matrix was used to perform the 
integration over the element volume (see Section 5.1.2). Thus the element stiffness matrix 
w r as obtained as the averaged stiffness of the different materials in the element. Periodic 
displacement and traction boundary conditions were imposed between corresponding nodes 
on opposite faces of the unit-cell. The numerical procedure to compute the element stiffness 
matrix and to impose periodic BC’s are discussed in Chapter 5. The skyline solver (Bathe, 
1982) was used to solve for the nodal displacements. The computed elastic constants for the 
seven examples are listed in Tables 3-5. 
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The fiTE %- 10 code was first checked by computing the elastic constants for an 
isotropic medium (Table 3). Then the code was implemented for a bimaterial medium. The 
bimaterial medium consisted of two different layers of equal thickness in the xy-plane 
alternating in the z-direction (Fig. 7). The effective Young’s moduli, Poisson ratios and 
CTE’s of the bimaterial medium were derived exactly, as described below. The constitutive 
relation (considering only the normal stresses) for each layer is defined as: 
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where the superscript refers to the layer number. To derive Cjj, Cj 2 and C;j for the bimaterial 
medium, e xx M was assumed as one; and e yy M and £ ZZ M were assumed to be zero. The 
assumption of £ zz M -0 and the fact that the layers are of equal thickness imply that e u 1 = 
-Szz 2 . The following constraints were, in addition, imposed across the bimaterial interface 
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From the above interfacial constraints and Eqn. (10), the stresses in each layer were 
computed. The stresses in the layers were volume-averaged to yield the corresponding 
macrostresses, i.e., o^, o yy M and o zz M . Since Exx M was equal to unity, the computed 
macrostresses were identical to the stiffness coefficients in the first column, namely, Cj], C ]2 
and C] 3 . A similar procedure was followed to find the remaining stiffness coefficients and 
CTE’s for the bimaterial medium. The inplane shear modulus of the bimaterial medium was 
computed as G xy = (G/ + G2)/2 knowing that the shear strain was uniform in both layers. The 

2GjG 2 

isostress assumption was used to derive the transverse shear modulus as G xz = ^ 
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It was found that^TEft-10 results were identical to the elasticity results for the bimaterial 
medium (Table 3). 


Table 3. Continuum properties for examples 1 and 2 using finite elements. 




E x , Ey 
(GPa) 

Ez 

(GPa) 

G x z> 
Gy z 

(GPa) 

Gxy 

(GPa) 

v xz> 

Vyz 


a x c , a y c 

X 

lO^C 

a-z c 

X 

io-Vc 

Example 1 
(isotropic 

fiTE%- 10 
(FEA) 

10 

10 

3.85 

3.85 

0.300 

0.300 

10 

10 

medium) 

exact 

solution 

10 

10 

3.85 

3.85 

0.300 

0.300 

10 

10 

Example 2 
(bimaterial 

fiTEx- 10 
(FEA) 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 

medium) 

exact 

solution 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 


Table 4. Continuum properties for examples 3 and 4 using finite elements. 




El 

(GPa) 

Ej 

(GPa) 

ESS 

Gjt 

(GPa) 

vlt 

VJT 

'Jj u 
a 2 

X 

dj 

xlO -6 /°C 

Example 3 

/uTEx-10 

(FEA) 

63.55 

36.48 

12.93 

9.94 

0.300 

0.232 

15.74 

40.79 

(uni 

directional 

composite) 

rule of 
mixt./ 
Halpin- 
Tsai eqns. 

64 

34.55 

11.26 

13.29 

0.300 

0.300 

15.63 

55.11 

Example 4 

fiTEy-lO 

(FEA) 

43.12 

18.15 

5.59 

3.92 

0.242 

0.222 

7.40 

25.44 

(uni 

directional 

composite) 

rule of 
mixt./ 
Halpin- 
Tsai eqns. 

43.40 

14.79 

4.45 

5.91 

0.260 

0.252 

6.77 

34.24 


Table 4 presents the elastic constants and CTE’s for the two unidirectional composite 
examples. The unidirectional composite unit-cell is shown in Fig. 8. The unidirectional 
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composite properties were compared with available analytical solutions. The rule of 
mixtures formulae were used to predict El and Vlt ; the Halpin-Tsai equations (Halpin and 
Tsai, 1969) for Ej, Glt and vjj and Schapery’s expressions (Agarwal and Broutman, 1990) 
for the thermal coefficients ai and aj. To compare Gjt, transverse isotropy was assumed in 

E t 

the T-T plane and computed as n - r. However, the finite element results show that 

Z(i + Vjj) 

the assumption of isotropy in the transverse plane is not valid. The assumption would give 
Gtt as 14.81 GPa and 7.43 GPa for examples 3 and 4 respectively, thus grossly 
over-estimating the finite element results. The relations for El, vlt and ai are exact when 
the poisson ratios are identical for the fiber and the matrix. 


Table 5. 


Continuum 


properties for examples 5, 6 and 7 using finite elements. 




Ex, Ey 
(GPa) 

Ez 

(GPa) 

Gxz j 
Gy Z 

(GPa) 

Gxy 

(GPa) 

V xz, 

v yz 

Vxy 

a x c ,a y c 

xlO" 6 / 0 

C 

a z c \ 

io-Vc 

Example 5 

(plain- 

weave) 

fiTEx- 10 
(FEA) 

11.81 

6.14 

1.84 

2.15 

0.408 

0.181 

28.36 

79.57 

Dasgupta 

results 

14.38 

6.25 

1.94 

3.94 

0.463 

0.167 

22.50 

86.00 

Example 6 

(plain- 

weave) 

/iTEx- 10 
(FEA) 

53.61 

10.88 

4.41 

4.72 

0.365 

0.128 

1.56 

22.71 

TEXCAD 

64.38 

11.49 

5.64 

4.87 

0.396 

0.027 

1.33 

20.71 

Example 7 

(5-hamess 

weave) 

fiTEx- 10 
(FEA) 

64.51 

11.33 

4.45 

4.85 

0.329 

0.047 

1.55 

22.03 

TEXCAD 

66.33 

11.51 

4.93 

4.89 

0.342 

0.034 

1.46 

21.24 


Figures 9 and 10 illustrate the weave patterns for the plain-weave (examples 5 and 
6) and the 5-hamess satin weave (example 7) respectively. The properties for example 5 were 
compared with Dasgupta’s results for an overall fiber volume fraction ( Vf ) of 0.26. The yam 
properties were not specified in Dasgupta et al. (1990). So the rule of mixtures and 
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Schapery’s expressions were used to obtain the yam properties from the given fiber and 
matrix properties. The^IE^-l 0 results for examples 6 and 7 were compared with TEXCAD 
—an approximate analytical method developed by Naik (1994). In both the examples the 
overall fiber volume fraction was 0.64. The unit-cell was divided into 13x13x7 elements for 
all three examples. The elastic constants for the three textile composite examples are 
presented in Table 5. It must be noted that iiTE %- 10 will marginally under-predict the 
stiffness moduli — since the yam cross-section in the numerical model is approximated as a 
polygon inscribed within the actual cross-sectional area. Consequently the yarn/fiber 
volume fraction in the numerical model will always be less than the theoretical volume 
fraction. 


2.2 Stress Gradient Effects 

The methods explained in Section 2. 1 assume that the unit-cells exist in all the three 
directions. This will be true in the case of thick textile composites. However there are many 
applications in which thin composites are used. In fact in order to take advantage of the 
properties of composites, the structures have to be made of thin plate like members with 
stiffeners for load transfer. In such cases there will be fewer unit-cells in the thickness 
direction. Thus the free surface effects will be predominant. There will be severe stress 
gradients through the thickness, and they will have an influence on the apparent stiffness and 
strength of the structure (Marrey and Sankar, 1993a; Sankar and Marrey; 1993a). 

The following simple example will illustrate the stress gradient effects on stiffness. 
Consider a layered medium consisting of alternating layers of materials of equal thickness 
with Young’s moduli Ej and E 2 respectively (Fig. 12a). Any micromechanical model would 
predict that the medium can be considered as a homogeneous orthotropic material at 
macroscale and also the effective Young’s modulus in the longitudinal direction is 
(E] +E2)/2, and there would not be any bending-stretching coupling in the principal material 
direction. However, if we consider a bimaterial beam consisting of the same two materials 
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(Fig. 12b), we will find that there is a bending-stretching coupling, and also the flexural 
rigidity cannot be predicted from the Young’s modulus of the homogeneous orthotropic 
medium and the total beam thickness. The bimaterial beam has properties and behavior 
different from the corresponding infinite medium. This phenomenon is observed in the 
transverse shear behavior also (Sankar and Marrey, 1993b). A similar behavior is also 
expected in thin textile composites where there are fewer unit-cells in the thickness 
directions, and the unit-cells are not subjected to a macroscopically homogeneous state of 
deformation as assumed in Section 2.1. 



Figure 12. Example to explain stress gradient effects. 

(a) Layered medium; (b) Bimaterial beam. 


One method of overcoming the above difficulties in thin textile composites is to 
model the composite as a plate/beam, and compute the structural stiffness properties (e.g., 
[A], [ B ] and [ D ] of the plate) directly from the unit-cell analysis instead of the continuum 
stiffness properties such as Young’s modulus, shear modulus etc. In the following sections 
we illustrate these concepts — first for a thin textile composite modeled as a beam (Sankar 
and Marrey, 1993b; Marrey and Sankar, 1993b) and then for a textile composite plate. The 
purpose of the beam model is to present the issues involved in computing the structural 
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stiffness coefficients as opposed to continuum elastic constants. Further the periodic BC’s 
are different from those in the continuum model. 


2.3 Unit-Cell Analysis for Beam Thermo-Mechanical Coefficients 


In this section, a procedure for finding the equivalent flexural stiffness properties and 
thermal coefficients of a textile structural composite beam is described. The textile 
composite beam is assumed to be in the xz-plane with unit-cells repeating in the x-direction. 
A state of plane strain parallel to the xz-plane is assumed. On the macroscale it is assumed 
that the beam is homogeneous and its behavior can be characterized by the following beam 
constitutive relation : 
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where P, M and V are the axial force, bending moment and transverse shear force resultants, 
respectively; [X] is the symmetric matrix of beam stiffness coefficients; e q m , x M and y 0 M are 
the midplane axial strain, curvature and transverse shear strain respectively; ct/>, and ay, 
respectively, are the beam thermal expansion, thermal bending and thermal shear 
coefficients. The midplane deformations are related to the midplane axial displacement 
transverse displacement w, and rotation ip as: 
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( 13 ) 


Actually,#/;, #72, #22 and K33 are similar to the laminate stiffness coefficients A//, Bjj, Djj 
and x 2 A55 f respectively. There is no equivalence for Kjj and K23 in the laminate theory, 
because the layers are assumed to be orthotropic, and they are rotated about the z-axis only. 
However, such a coupling between inplane deformations and transverse shear deformation 
may exist in textile composites as the fibers are inclined to the xy-plane unlike in the 
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laminates. The beam constitutive relation in Eqn. (12) can also be expressed in terms of 
compliance coefficients as: 
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2.3.1 Beam Steady State Loading Conditions 



The continuum unit-cell analysis assumes that all the unit-cells are subjected to 
identical stress and strain fields, for a given state of loading. This is true in the case of 
constant axial force (P) and constant bending moment (. M ) in the beam (Figs. 13a and 13b). 
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However, when a shear force (V) is applied to the beam, the shear force will give rise to 
building up of bending moment at every cross-section, such that V = -( dM/dx), This bending 
moment varies linearly over the length of the beam violating the assumption of 
homogeneous deformation. A state where the unit-cells are subjected to identical 
deformation under a shear force can be created by adding a couple periodically (Fig. 13c), 
or by having shear tractions on the top and bottom surfaces to cancel the bending moment 
continuously (Fig. 13d). In both cases traction-free conditions are violated on the top and 
bottom surfaces of the beam. As will be seen later, this situation creates difficulties in 
estimating the shear stiffness of the beam accurately. 

2.3.2 Unit-Cell Boundary Conditions 

The left end of the unit cell, x=0 (Fig. 14) are subject to minimum support constraints 
to prevent rigid body translation and rotation. The top and bottom surfaces of the beam are 
assumed to be traction free. The edges x=0 and x~L have identical nodes in the finite element 
model and periodic BC’s are enforced between these nodes. Three linearly independent 
deformations were applied to the unit-cell, in order to find the beam stiffness matrix [ K] 
namely, 

Case (i) unit axial strain ( £% = l,x M = 0 ,y% — 0) 

Case (ii) unit curvature accompanied by transverse deflection such that the transverse 
shear strain was zero = 0,x M = 1 = 0) 

Case (iii) unit transverse shear strain ( e % = 0,x M = 0 = 1) 

The periodic displacement constraints applied for each deformation case are presented in 


Table 6. 
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Table 6. Periodic displacement BC’s for beam unit-cell. 



u(L,z)-u(0,z) 

w(L,z)-w(0,z) 

mm 

Case i. unit axial strain ( £o M =l ) 

L 

0 

0 

Case ii. unit curvature (x M =l) 

Lz 

-L 2 /2 

0 

Case iii. unit shear strain (yo M =l) 

0 

L 

0 


Consider the second deformation case, where a curvature is applied to the unit-cell. Let ’m’ 
and ’n’ be a set of corresponding nodes on the left and right ends of the unit-cell. The periodic 
BC’s imposed between the nodes ’m’ and ’n’ are shown in Fig. 14. Applying a curvature to 
the unit-cell also induces transverse shear strain due to the term in the expression for y q 
(Eqn. 13). The difference in displacement in the z-direction is therefore imposed to make the 
macroscopic transverse shear strain equal to zero. 
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2.3.3 Determination of Beam Stiffness Coefficients 


The three linearly independent deformations described in the preceding section are 
applied to the unit-cell. The temperature difference (AT) is set to zero for all three 
deformations. For each case, the axial force P, the bending moment M c and the shear force 
resultant V are computed from the nodal forces at the ends of the unit-cell. As explained 
earlier, the axial force and shear force are constant through the unit-cell, but the bending 
moment varies linearly across the unit-cell length. Therefore the bending moment at the 
center of the unit-cell is used as the reference moment, which is obtained by averaging the 
bending moment resultants at the ends of the unit-cell. From the force resultants, a 
pseudo-stiffness matrix [fc] can be computed that relates the forces and deformations as 
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For example, kjj, JC 21 and kjj will correspond to the values of P, M c and V obtained for the 
case of unit axial extension of the unit-cell (Case i). The pseudo stiffness matrix, in general 
will not be symmetric since it does not relate corresponding forces and deformations 
(conjugate quantities, product of which yields an energy term) — it can be rather considered 
as a matrix of influence coefficients. The inverse of [&], denoted by [ s ] has some significance. 
The matrix [s] is defined as 
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(16) 


We know that V=0 is a steady state loading condition such that M C =M, where M is the 
constant bending moment along the unit-cell. Noting that only the last column of [s] 
multiplies with V, and comparing Eqns. (14) and (16), we can conclude that the first two 
columns of [S] and [ 5 ] must be identical to each other. Since [S] is symmetric (Sj 3 - 831 , 
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S23=S32 h we can determine all but S 33 of the beam compliance matrix [S]. An alternative 
approach to predict the shear compliance S 33 is explained in the next section. 

2.3.4 Determination of Beam Shear Stiffness 

The difficulty in estimating S 33 (or K 33 ) is associated with the inability to create a 
state of deformation such that only V is present. The shear modulus of the textile composite 
G xz can be computed by assuming that the unit-cells span the material in the z-direction also. 
In that case periodic BC’s are imposed between the top and bottom surfaces of the unit-cell, 
and the unit-cell is deformed transversely. Then, there will be shear tractions on the top and 
bottom surfaces of the beam — in fact this situation would correspond to Fig. 13(d). One may 
surmise that a shear correction factor x 2 could be found such that K 33 ~K 2 G xz h. But a simple 
bimaterial beam example will show that the shear stiffness can be grossly underestimated. 

Consider a bimaterial beam with layers of equal thickness (h/2). By performing a 
continuum unit-cell analysis the shear modulus of the material is found to be equal to 2Gj G 2 
/ (Gj+ G 2 ). Assuming G]/G 2 - 10, G xz - 0.182 Gj. However the actual shear stiffness K 33 
is equal to x 2 (Gj+ G 2 ) h/2, which means that the apparent shear stiffness G xz = (Gj+ G 2 )/? 
= 0.55G 1 . This is about three times the previous estimate. This discrepancy is due to differing 
assumptions regarding the constancy of shear stress or shear strain. The continuum unit-cell 
analysis imposes constant shear stress in the two materials, and hence the shear compliance 
of the composite is the average of the compliances of the constituent materials. This is true 
when there are a large number of unit-cells in the z-direction. In a bimaterial beam, however, 
the shear strain is almost constant in the two layers, and hence the shear stiffness is the 
average of the shear stiffness of the individual layers — which is consistent with the method 
of computing A 55 in the lamination theory. This illustrates the need for special procedures 
for predicting the shear stiffness of thin textile composite beams. 

To overcome this problem we use an energy approach to compute the shear stiffness 
of a thin textile composite. We model a beam of length 2L consisting of two unit-cells, with 
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plane strain finite elements. The beam is subject to boundary conditions corresponding to 
pure shear strain (third boundary condition in Table 6 ). The top and bottom surfaces of the 
beam are traction free. The shear strain energy over a length L in the middle of the beam, U s , 
is computed from the finite element results as 



where x^ xz and y^ xz are the shear stress and shear strain at the center of the i th element and 
Ai is the area of the i th element. The summation is performed over all elements located in a 
length L in the middle of the beam. Next, the shear strain energy over the same length L is 
computed using the beam formula: 

U s = ( r)vy 0 i = (Ijw + S 23 M c + S 33 V)L (18) 

In the above equation P, M c and V can be obtained from the finite element results. The 
coefficients S 73 and S 23 have already been estimated. Then by equating the shear strain 
energy quantities in Eqns. (17) and ( 1 8), the only unknown, that is S 33 , can be evaluated. The 
choice of two unit-cells to perform the above analysis deserves an explanation. When this 
was tried with one unit-cell for the cases of isotropic and bimaterial beam, the results were 
not good. It was mentioned earlier, that the application of a shear force would result in a 
couple at the unit-cell ends to continuously cancel the effect of the bending moment created 
along the beam length. This couple is manifested in the finite element results as concentrated 
forces at the four corners of the unit-cell, thus creating severe stress concentrations. When 
two unit-cells are used in the model, the stress concentrations remain in the corners of the 
beam, but their effects diminish in the middle portion of the beam. As will be seen in the later, 
the two unit-cell method gave very good K 33 for both isotropic and bimaterial beams. 

2.3.5 Determination of Beam CTE’s 

The procedure for determining the beam thermal expansion coefficients (Marrey and 
Sankar, 1993b) is as follows. The beam unit-cell is subject to a uniform temperature 
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difference given by AT = T 0 . The deformations in the beam are restrained by setting 
eo M =x M =yo M =0. Then the beam constitutive relation, Eqn. (12), will reduce to 
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The axial force P ; the bending moment resultant at the center of the unit cell M c and the shear 
force V are computed from the nodal forces at the ends of the unit-cell. Then the beam CTE’s 
can be estimated from the expression. 
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2,3.6 Results and Discussion 


The procedures described above were implemented for the following cases: 

(a) an isotropic beam 

(b) a bimaterial beam with isotropic layers of equal thickness 

(c) a plain weave textile composite beam where the yam is assumed to be transversely 
isotropic and the matrix is isotropic. 


Table 7. Constituent material properties for beam examples. 


isotropic beam 

E- 10 GPa, v - 0.30, a - lOxlO" 6 /^ 

bimaterial beam 

E j = 70 GPa, vj = 0.33, aj = 23X10" 6 /°C 
E 2 = 3.5 GPa, v 2 = 0.35, a 2 = 60X10- 6 PC 

plain-weave 
textile beam 

yam: 

E; = 159 GPa, E 2 = 10.9 GPa, G n = 6.4 GPa, 
v ]2 = 0.38,v 2 j = 0.38, aj = 0.045x10^ PC, a 2 = 20.2X10" 6 PC 
where the yam direction is parallel to the i-axis and 23-plane 
is the plane of isotropy. 

matrix: 

En, = 3.5 GPa, v m = 0.35, a m = 60x10-6 /°C. 
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The properties of the constituent materials for all the cases are listed in Table 7. The 
dimensions of the unit-cell and the yarn architecture for the textile beam were taken from 
Yoshino and Ohtsuka (1982). The same unit-cell dimensions (length of 3.6 mm and height 
1.8 mm) were also used for the isotropic and bimaterial cases. The unit-cell of the beam was 
discretized using eight-node isoparametric plane strain finite elements. The finite element 
mesh for the isotropic unit-cell and the plain weave unit-cell were identical except that 
different material properties were used. 


Table 8. Comparison of beam stiffness coefficients and CTE’s (SI units). 




Kn 

&12 

K 22 

K33 

a P /°C 

a M /°C 

isotropic 

beam 

unit-cell 

analysis 

19.78x106 

0 

5.35 

5.96x106 

lOxlO^ 6 

0 

beam 

theory 

19.78x106 

0 

5.34 

5.77x106 

10x10-6 

0 

bimaterial 

beam 

unit-cell 

analysis 

74.29x106 

30.20x103 

20.06 

8.47x106 

30.73x1 CH 

-14.62x10-3 

beam 

theory 

74.29x106 

30.20x103 

20.06 

8.62x106 

30.74x10-6 

-14.63x10-3 

textile 

beam 

unit-cell 

analysis 

27.76x106 

0 

5.41 

9.21x106 

12.66x10-6 

-24.12x10-6 

mosaic 

model 

71.48x10-6 

0 

8.13 

8.14x106 

4.39x10-6 

0 


Note; K] 3 , K 23 and ay are zero for all cases 


The deformed unit-cells under various independent loading conditions are shown in 
Figs. 15 to 17. The stiffness and thermal coefficients for the three beams are shown in Table 
8. The results for the isotropic and bimaterial beams were compared to exact beam theory 
solutions. Exact shear correction factors — 0.833 for the isotropic beam and 0.555 for 
bimaterial beam (Whitney, 1 973) were used in the beam theory solution to compute the shear 
stiffness given in Table 8. It can be seen from Table 8 that the beam unit-cell analysis is able 
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to predict the axial and bending stiffness coefficients {Kjj and K 22 ) very accurately. As 
expected the shear stiffness ( K 33 or A 55 ) predictions have errors, but they are very minimal. 
It can be noticed that the comers of the unit-cell are severely deformed (Figs. 15c, 16c and 
1 7d), when the unit-cell is subject to constant shear strain leaving the top and bottom surfaces 
traction free. However when shear tractions are allowed on the top and bottom surfaces of 
the unit-cell, the distortions at the comers disappear (Figs. 15d, 16d and 17e). Then what is 
obtained is the shear modulus G xz and not the beam shear stiffness. The shear modulus of 
the plain weave beam was found to be 3 .07 GPa. This would yield the apparent shear stiffness 
as G xz h =5.53 x 10 6 Nm -1 ; whereas the actual shear stiffness is 9.2 1 x 10 6 Nm -1 ( K 33 in Table 
8 ). The Young’s modulus of the textile beam E x may be extracted from Kjj , as Ku/h, which 
would yield E x = 15.42 GPa. If this value of E x were used to predict the flexural stiffness of 
a homogeneous beam as Du = Exh 3 /12, we would obtain Djj as 7.50 Nm- — whereas the 
actual flexural stiffness is equal to 5.41 Nm. The same idea holds for the beam thermal 
coefficients also. The beam CTE’s ap, (Xm and ay cannot be predicted from the 
corresponding continuum CTE’s. Table 9 shows the discrepancy, for the plain weave 
example, between the beam CTE’s obtained directly and the beam CTE’s obtained from the 
corresponding continuum CTE’s. It may be noted that the continuum model would always 
predict the thermal expansion coefficient ap as a x , and the thermal bending coefficient a m 
as zero. This underscores the importance of the present analysis for predicting the beam 
stiffness properties for a thin textile composite directly. 


Table 9. Comparison of textile CTE’s. 



CTE’s from 

CTE’s from 

% error 


beam model 

continuum model 


ap x ICC 6 PC 

12.66 

11.30 

-10.46 

aM x 10 -6 /°C/m 

-24.12 

0 

00 

ay PC 

0 

0 

0 
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Figure 15. Deformed unit-cell of the isotropic beam, (not to scale) 

(a) unit extensional strain; (b) unit curvature; (c) unit shear strain, top and 
bottom surfaces are traction free; (d) unit shear strain, tractions allowed on 
top and bottom surfaces. 
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The textile beam stiffness coefficients were also estimated using a procedure similar 
to the mosaic model (Ishikawa and Chou, 1982a). They are compared with the coefficients 
obtained from the unit-cell analysis (Table 8). The idealization made in the mosaic model 
is shown in Figure 18. The unit cell is divided into five segments, with each segment modeled 
as a cross-ply laminate with a stacking sequence which best represented the yarn architecture 
within that segment. The stiffness matrix of each segment was computed using laminate 
analysis. Then the compliance of the textile beam was computed as the length-weighted 
average of the compliance of the five segments. 

5 

[5] = 4 ) 2 L k [S] k (21) 

7 k= 1 

From Table 8, it can be seen that the mosaic model predicts K 33 reasonably well. The reason 
for the lack of agreement in Kjj and K 22 can be attributed to the fact that a major portion of 
the yam is modeled as a 0° laminate in the mosaic model, which tends to over-predict the 
axial and flexural stiffnesses. 



Figure 18. Mosaic model of the textile unit-cell. 
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2.4 Plate Thermo-Mechanical Coefficients 


The previous section described the concepts involved in modeling a thin textile 
composite structure as a., homogeneous beam to predict the beam stiffness coefficients and 
CTE’s. The procedure is extended in this section, utilizing three-dimensional finite element 
analysis, to model the thin composite as a plate (Sankar and Marrey, 1992) to determine the 
corresponding thermo-elastic coefficients. 

The plate is assumed to be in the xy-plane with unit-cells repeating in the x - and 
y-directions. The lengths of the unit-cell in the x- and y - directions are assumed to be a and 
b respectively and the unit-cell thickness as h. On the macroscale the plate is assumed to be 
homogeneous and the plate behavior is characterized by the plate constitutive relation: 


Ny 

Nxy 

M x 

My 

M*y 



(22) 


where £io M , yio M and K[ M are the midplane axial strain, shear strain and curvature; aj? and 
Pi p are the plate thermal expansion and bending coefficients; fy and Af; are the axial force 
and bending moment resultants respectively in the homogeneous plate. The plate stiffness 
matrix comprises of the [A], [B] and [D] sub-matrices, which are the plate extensional 
stiffness, coupling stiffness and bending stiffness matrices respectively. The plate stiffness 
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matrix can be expanded as: 


A B 
B D 


A 

ll 

^12 

^16 

B 

ll 

*12 

*16 

A 

12 

^22 

^26 

B 

12 

*22 

*26 

A 

16 

^26 

^66 

B 

16 

*26 

*66 

B 

11 

*12 

*16 

D 

11 

**12 

0,6 

B 

12 

*22 

*26 

D 

12 

D 22 

^26 

B 

16 

*26 

*66 

D 

16 

D 26 

0 6 6 


(23) 


The midplane strains and curvatures are related to the midplane displacements and rotations 
as: 


M _ 
£ x0 ~ 


X 


M _ 


Bu a 
~Bx ’ 

11 

to 

dv 0 

dy' 

II 

Bu 0 dv 0 

By 3x 

Bx ’ 

II 

dlpy 

By ’ 

II 

Bxpx dtyy 

By Bx 


2.4.1 Unit-Cell Boundary Conditions 


(24) 

(25) 


The plate thermo-mechanical properties are obtained by modeling the unit-cell with 
eight-node brick elements and subjecting the unit-cell to six linearly independent 
deformations. The six deformations are given by: (1) unit e ^ maintaining the rest of the 
macroscopic strains and curvatures as zero; (2) unit e such that remaining strains and 
curvatures are zero; and similarly (3) unit yj^; (4) unit x^\ (5) unit x™\ (6) unit x%. In the 

last three cases (non-zero curvatures) the deformation was accompanied with a transverse 
deflection such that the transverse shear strain was zero (Table 10). 
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Table 10. Periodic displacement BC’s imposed on the lateral faces of the plate unit-cell. 



The unit-cell is subject to minimum support constraints to prevent rigid body rotation 
and translation (Fig. 19). The top and bottom surfaces of the plate are assumed to be free of 
tractions. The faces x—0 and x=a have identical nodes in the finite element model, and so 
do the pair of faces y=0 and y=b. The identical nodes on opposite faces of the unit-cell are 
constrained to enforce the periodic BC’s. The traction boundary conditions on the lateral 
faces of the unit-cell are: 
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Fi(a,y,z) = - Fj(0,y,z) 

F £x, b,z)= - Ffx, 0, z) i = x ’ y ' z (26) 

The periodic displacement BC’s enforced for each mode of deformation are presented in 
Table 10. 


2.4.2 Determination of Plate Stiffness Coefficients and CTE’s 


The six linearly independent deformations are applied to the unit-cell such that only 
one of the macroscopic strains or curvatures is non-zero (first six cases in Table 10). The 
temperature difference is set to zero for all six cases. It must be noted that the applied 
deformations must ensure that the transverse shear strains, y xz M and y yz M are zero where 


v* = ** + f 


(27) 


The force and moment resultants can be obtained by one of the following two methods. In 
the first method, the resultants are computed by averaging the nodal forces on each face of 
the unit-cell. For example, on the face x~a the force and moment resultants are computed 
using the relations: 



where F x ^ and F y b) are the nodal forces in the x and y directions at the i th node and V is 
the total number of nodes on the face. The force and moment resultants can also be computed 
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by averaging the microstresses over the unit-cell volume. Then the resultants on the face jc -a 
are obtained as 


c (*,y, z)dV 


p 

_L 

ab 


Txy(x,y,z)dV 


M x ” ^ J zo xx {x,y,z)dV Mxy = ~ J zr xy (x i y,z)dV (30) 

V V 

Substituting the values of the deformation and the force resultants in the plate constitutive 
relation, Eqn. (22), the stiffness coefficients in the column corresponding to the non-zero 
deformation can be computed. This procedure is repeated for other deformation components 
to obtain all the stiffness coefficients. 

To predict the CTE’s, the plate unit-cell is subject to a uniform temperature 
difference, given by A T — Tq. In the finite element model, periodic displacement BC’s are 
applied such that all six components of the deformation are zero (seventh case in Table 10). 
The averaged force and moment resultants are computed using one of the procedures 
described above. The thermal expansion coefficients aP, and thermal bending coefficients 
pp are then obtained from the relation: 



2.4.3 Results for Plate Stiffness Coefficients 

The plate [A], [5], [D] matrices and CTE’s were found for the seven examples listed 
in Table 2 by implementing the finite element cod q^iTEx-IO. As a preliminary check, the 
code was executed for an isotropic plate and compared with the plate properties using 

Eh 

lamination theory (for one-ply). For example, Ajj was calculated as r- and Djj as 

(1 - v z ) 

ril 3 

. Then the properties were computed for a bimaterial plate using the finite 
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element code. The plate properties for the isotropic and bimaterial cases are presented in 
Tables 1 1 and 12 respectively. The bimaterial plate properties were also computed using the 
lamination theory for two plies, and from the homogenous 3-D elastic constants computed 
in the previous section (Section 2.1). For example the coefficient Djj is obtained from the 

E M h 3 

3-D elastic constants as Di , = — r- . The finite element results for the bimaterial 

11 12(1 - v" 2 ) 

case were exact, i.e., identical to the results obtained with the two-ply lamination theory. The 
[ D ] matrix computed from the bimaterial 3-D constants was found to be in good agreement 
with the two-ply lamination theory only because both the layers were equal in thickness. This 
is a special case, and in general, the [D] matrix obtained from the 3-D elastic constants will 
be different from the two-ply lamination theory results. 

The plate properties for the unidirectional composite examples are presented in Table 
13 and for the textile examples in Table 14. In all the examples it was found that the plate 
properties, especially , {a p } and {fi p } could not be predicted from the corresponding 

3-D elastic constants. 


Table 11. Non-zero [A], [B] and [D] coefficients for example 1 (isotropic plate) using 
finite elements 



An 

xlO 6 

An 

xlO 6 

A 22 

xlO 6 

A 56 

xlO 6 

D 11 

xlO -3 

D]2 

xlO -3 

D 22 

xlO -3 

A 56 
xitr 3 

a x P, a y P 
xl0^/°C 

HTE X - 10 
(FEA) 

2.810 

0.843 

2.810 

0.983 

15.320 

4.606 

15.320 

5.358 

10 

lamination 

theory 

2.810 

0.843 

2.810 

0.983 

15.310 

4.593 

15.310 

5.358 

10 


Note: [A], [£] and [/>] coefficients in SI units 
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Table 12. Non-zero [A], [ B ] and [D] coefficients for example 2 (bimaterial plate) using 
finite elements 



All, A 22 
xlO 6 

A, 2 
xlO 6 

A> 6 
XlO 6 

Bn, B 22 
xlO 3 

Bl 2 
xlO 3 

^66 

xlO -3 

fiTEx- 10 
(FEA) 

9.832 

2.043 

3.895 

-0.563 

-0.108 

- 0.228 

lamination theory 
for two plies 

9.832 

2.043 

3.895 

-0.563 

-0.108 

- 0.228 

lamination theory 
using 3-D elastic 
constants 

9.844 

2.048 

3.899 

0 

0 

0 



Du , D 22 
xlO- 3 

B>12 

XlO -3 

D 66 

xlO -3 

a x p , a y P 
xlO _6 /°C 

/°C/m 

juTEx - 10 
(FEA) 

53.590 

11.149 

21.220 

17.800 

0.170 

lamination theory 
for two plies 

53.573 

11.131 

21.220 

17.814 

0.170 

lamination theory 
using 3-D elastic 
constants 

53.762 

11.183 

21.293 

8.190 

0 


Note: [A], [5] and [ D ] coefficients in SI units 
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Table 13. Non-zero [A], [B] and [D] coefficients for examples 3 and 4 (unidirectional 
composite) using finite elements 



Mi 

An 

A22 

A66 

xlO 6 

xlO 6 

xlO 6 

xlO 6 


Example 3 

yUlEft-lO 

(FEA) 

0.690 

0.149 

0.496 

0.177 


Halpin-Tsai Eqns. 
and lamination 
theory 

0.673 

0.109 

0.363 

0.113 

Example 4 

[iTEx-\0 

(FEA) 

0.452 

0.062 

0.285 

0.114 


Halpin-Tsai Eqns. 
and lamination 
theory 

0.444 

0.039 

0.151 

0.045 



Dn D]2 D22 D $6 a x Px dyPx 

xlO -6 xlO” 6 xlO" 6 xlO- 6 10- 6 10- 6 


Example 3 

tiTEx - 10 
(FEA) 

3.589 

0.596 

1.980 

0.947 

15.489 

26.184 


Halpin-Tsai 
Eqns. and 
lamination 
theory 

5.606 

0.908 

3.026 

0.939 

15.625 

55.112 

Example 4 

fiTEx - 10 
(FEA) 

2.256 

0.224 

0.873 

0.568 

7.378 

13.188 


Halpin-Tsai 
Eqns. and 
lamination 
theory 

3.702 

0.328 

1.262 

0.371 

6.774 

34.239 


Note: [A], [5] and [D] coefficients in SI units 
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Table 14. Non-zero [A], [ B ] and [D] coefficients for examples 5, 6 and 7 using finite 
elements 



A n , A 2 2 A] 2 A6 6 Bjj xio 3 

xlO 6 xlO 6 xlO 6 


Example 5 

07^-10 

(FEA) 

2.681 

0.565 

0.489 

0 


lamination theory 
using 3-D constants 

2.783 

0.503 

0.490 

0 

Example 6 

fiTE %- 10 
(FEA) 

12.090 

3.470 

1.223 

0 


lamination theory 
using 3-D constants 

13.938 

1.787 

1.208 

0 

Example 7 

JuTE X - 10 
(FEA) 

14.683 

1.351 

1.210 

0.495* 


lamination theory 
using 3-D constants 

16.531 

0.770 

1.239 

0 



Dll, D 22 &12 

xlO -3 xlO“ 3 



Ox?, CLyP p x P 
xlO" 6 /°C /°C/m 


Example 5 

fiTEx - 10 
(FEA) 

5.687 

1.518 

1.577 

27.465 

0 


lamination 
theory using 
3-D constants 

12.054 

2.177 

2.124 

28.363 

0 

Example 6 

juTEx - 10 
(FEA) 

41.695 

0.373 

5.879 

1.480 

0 


lamination 
theory using 
3-D constants 

75.942 

9.734 

6.582 

1.556 

0 

Example 7 

pTEx - 10 
(FEA) 

90.072 

1.123 

6.149 

2.910 

-0.037* 


lamination 
theory using 
3-D constants 

87.283 

4.195 

6.753 

1.550 

0 


* In example 7, # 22 = -Bn and p y P = -p x P 

Note: [A], [B] and [ D ] coefficients in SI units 










































































CHAPTER 3 

ANALYTICAL MODELS FOR THERMO-ELASTIC CONSTANTS 


The complex yarn architectures in a textile composite make numerical modeling of 
the unit-cell extremely difficult. Besides, the computational memory and run-time 
requirements for a detailed finite element analysis are enormous. As explained in the 
introductory chapter, there are several parameters that can be changed to alter the effective 
composite properties. These parameters may be the fiber material in the yam, fiber volume 
fraction in the yam (also called yam packing density), overall fiber volume fraction, preform 
architecture or the matrix material properties. This emphasizes the need for simple analysis 
procedures to predict the trend in variation of composite properties when one of the 
parameters is changed. These procedures will be of use to a designer in determining the 
optimum parameters for a certain application. 

Analytical methods are approximate because they assume certain forms for the state 
of stress and strain in the unit-cell. Averaging the stiffnesses or compliances of the matrix 
and the inclusion has long been used to estimate the bounds of effective elastic properties 
of the composite. Essentially the stiffness averaging assumes a state of uniform strain in the 
composite (isostrain), and compliance averaging assumes a state of uniform stress (isostress) 
in the matrix and inclusion. In fact the rule of mixtures expressions for estimating the 
effective properties of a unidirectional composite is based on such averaging schemes. Naik 
(1994) proposed an analytical method (TEXCAD) in which the yams are discretized into 
segments. Knowing the direction of the yam in each segment, the segment stiffness are 
computed using appropriate transformations. Then assuming a state of isostrain, the 
composite stiffness is obtained by volume-averaging the yam-segment stiffness and matrix 
stiffness in the unit-cell. This method seems to work, when there is multi-directional 
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reinforcement in the composite. However the method fails for composites with preferential 
yam reinforcement. For example, the transverse modulus and the inplane shear modulus for 
a unidirectional composite cannot be estimated using the above method. Another analytical 
model called the mosaic model, was proposed by Ishikawa and Chou (1982a). As discussed 
earlier, in the mosaic model the yam architecture is simplified to that of a cross-ply laminate 
and the lamination theory is used to predict the composite properties. 

The state of stress/strain in a textile composite, subjected to a uniform macrostress, 
is much more complex than that assumed in the above mentioned methods. We propose a 
scheme of selective averaging — called the selective averaging method (SAM) — in which 
both stiffness and compliance coefficients can be averaged selectively depending on a more 
realistic assumption of either isostress or isostrain. 

3.1 Selective Averaging Method fSAMl for Continuum Properties 

In this section, we describe SAM for estimating the effective elastic constants and 
CTE’s for a textile composite. Consider a rectangular hexahedron of dimensions axbxc 
as the unit-cell. The unit-cell is discretized into slices on the mesoscale, and elements on the 
microscale as shown in Fig. 20. To distinguish between the macrolevel, mesolevel and 
microlevel properties in this section, an over-tilde is used to denote the mesolevel properties, 
and a superscript ”M” is used to denote the macrolevel properties. For example, [C M ], [C] 
and [C] will represent the macrolevel, mesolevel and microlevel stiffnesses respectively. 

The objective here is to determine the coefficients of the effective stiffness matrix 
[C M ] as defined in Eqn. (2). To find the first column of the effective stiffness matrix, the 
unit-cell is divided into slices (mesolevel) of thickness dx parallel to theyz-plane (Fig. 20a). 
Each slice is further sub-divided into elements (microlevel) as shown in Figs. 20(b) and 
20(c). The unit-cell is subjected to a deformation such that all macrostrains except £xx M are 
equal to zero and £xx M =1- It is assumed that the mesolevel and microlevel strains, 
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corresponding to the zero macrostrains, are negligible. In other words, 

Sf - E t = £, = 0 i 1 (32) 



Figure 20. Hierarchy of discretization for a unit-cell to determine first and fifth 
columns of effective stiffness matrix. 

(a) unit-cell; (b) slice; (c) element. 


Assuming a state of isostrain within the slice (£**(*, y,z) = £(*) ) the average stiffness of 
a slice can be obtained as: 


C n (x) 


be 




C u (x,y,z) dy dz 


(33) 


where Cn(x,y,z) is the element stiffness coefficient referred to the unit-cell coordinates. The 
stiffnesses of the slices are averaged on the macrolevel based on the isostress assumption, 
i.e., Ojocix) = (jM . Then the first column of the effective stiffness matrix can be computed 
using the following two relations: 
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a 



( 34 ) 


r^M — 1 

11 abc 


c b a 

f f \ 4 &- 

J J J C n (x) 


C iX (x,y,z) dx dy dz (i = 2,..., 6) 


(35) 


z=0>’ = 0x=0 


A similar procedure can be implemented to determine the second and third columns of the 
macroscale stiffness matrix \C M ]. 

A slightly different averaging scheme is used when the unit-cell is subjected to shear 
strains on the macrolevel. Consider the case where the unit-cell is subjected to unit y yz at 
macroscale. The unit-cell is again discretized into slices and elements as shown in Fig. 20. 
It is assumed that all the other components of strain at the macrolevel, mesolevel and 
microlevels are zero. This can be expressed as: 


ef — = £. = 0 i * 4 (36) 

where £4 = y yz . We also assume that the shear stress is constant in a slice such that 
r yz {x, y, z) = T yz (x). The shear compliance of a slice can then be obtained by averaging the 
shear compliances of all the elements in the slice as: 

c b 

1 = j_ 

C M (x) bc 

The fourth column of the stiffness matrix Q 4 M is obtained under the assumption that all the 
slices are under a state of constant shear strain: 


J 

= 0 


1 


C44(*,y,z) 


dy dz 


(37) 


c b a 


c % = c^ixy^z) Ci4 ^ x,y ’ z ) ^ dy dz V = 


z = 0y = 0 x=0 


A similar procedure is used to determine the fifth and sixth columns of [ C M ]. 



52 


To determine the macroscale CTE’s, a uniform temperature difference {AT) is 
applied throughout the unit-cell. The unit-cell is constrained from expanding such that all 
the strain components on the macrolevel are zero. A state of isostrain is assumed in the 
unit-cell, implying that the mechanical strain components on the mesolevel and microlevel 
are also zero. This can be expressed as: 

ef = £,. = £, = 0 i = (39) 

Then the thermal constitutive relations on the macrolevel and microlevel will reduce to: 

\o M ) = - [C M ]{a c )AT 

1 J (40) 

{or} = - [C]{a}AT K ; 


The macrostresses may be computed by volume-averaging the corresponding microstress 
component as shown below: 

c b a 

[o M ] = ^ I I I {a} dx dy dz (41) 

z=0y==0x=0 

Then from Eqns. (40) and (41), we can compute the macroscale CTE’s as: 

f« c > =Az [CV {/} (42) 


where {/} is given by the expression: 


{/} 


c b a 



[C){a } dx dy dz 


z = 0>- = 0jc=0 


(43) 


3.2 Continuum Results using SAM 

A code called fxTE%-20 (pronounced as microtech) was written in FORTRAN 77 to 
implement SAM. The code was executed to estimate the thermo-elastic constants for the 
seven examples, whose constituent material properties are listed in Table 1 . Input to the code 
were the unit-cell dimensions, yam geometry information, constituent material properties, 
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and the number of divisions required to discretize the unit-cell in the y and z directions. 
The element stiffness matrix [C] was determined by computing the elasticity matrix for the 
material point at the geometric center of the element, and transforming it to the unit-cell 
coordinate system. The predicted macroscale stiffness matrix [ C M ], and consequently the 
macroscale compliance matrix will not be symmetric due to the approximate nature of the 
analysis. Therefore the macroscale compliance matrix was made symmetric by averaging 
the off-diagonal compliance coefficients. The macroscale elastic constants were computed 
by comparing the symmetrized compliance coefficients with that of a homogenous, 
orthotropic medium. The user manual for [iTE %- 20 is provided in Appendix- A. 

The results for a bimaterial medium (example 2) are given in Table 15. The results 
for example 1 are not listed, since it is obvious that SAM would predict the elastic constants 
for an isotropic medium exactly. The bimaterial medium consisted of two different layers 
of isotropic materials of equal thickness altematingly stacked in the z-direction. The elastic 
constants for the bimaterial medium were compared with an exact solution (derivation 
explained in Section 2. 1 .3), and with the previously computed finite element results. It can 
be observed that SAM marginally under-predicts the longitudinal and transverse Young’s 
moduli, while the inplane and transverse shear moduli are exact. 


Table 15. Continuum properties for example 2 using SAM. 




Ex, Ey 

(GPa) 

Ez 

(GPa) 

G xz , 

Gyz 

(GPa) 

Gxy 

(GPa) 

Vxz, 

v yz 


a x c , a y c 

X 

lo-fyoc 

a z c 

X 

io- 6 /^ 

Example 2 
(bimaterial 
medium) 

[iTEx - 20 
(SAM) 

36.02 

8.72 

2.48 

15.23 

0.599 

0.183 

3.88 

52.20 

HTE%- 10 
(FEA) 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 

exact 

solution 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 
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Table 16 presents the SAM results for two cases of unidirectional composite 
(examples 3 and 4). The fiber and matrix had identical Poisson’s ratio in example 3, and 
different Poisson’s ratio in example 4. The SAM results were compared with the finite 
element results from the previous chapter, and with analytical solutions for unidirectional 
composite properties. The analytical expressions used were the rule of mixtures formulae 
for Ei and vu and the Halpin-Tsai equations (Halpin and Tsai, 1969) for Ej, Gu and Vjj- 
All of the unidirectional composite thermo-elastic constants but for Ej and aj were found 
to match well with the compared data. Table 17 compares the SAM results for three textile 
composites (examples 5, 6 and 7) with available results. In all three cases the thermo-elastic 
constants obtained by implementing SAM were in good agreement with the available results. 


Table 16. Continuum properties for examples 3 and 4 using SAM. 




El 

(GPa) 

Ej 

(GPa) 

Glt 

(GPa) 

Grr 

(GPa) 

V LT 

Vtt 

X 

to- 6 /°C 

aj 

x\Qr (> /°C 

Example 3 

(unidirect. 

composite) 

fiTE%- 20 

(SAM) 

64 

50.24 

10.45 

8.36 

0.341 

0.300 

12.41 

21.51 

fiTEx- 10 
(FEA) 

63.55 

36.48 

12.93 

9.94 

0.300 

0.232 

15.74 

40.79 

rule of 
mixt./ 
Halpin- 
Tsai eqns. 

64 

34.55 

11.26 

13.29 

0.300 

0.300 

15.63 

55.11 

Example 4 

(unidirect. 

composite) 

fiTEx- 20 
(SAM) 

43.35 

32.47 

4.13 

3.04 

0.245 

0.218 

7.40 

11.60 

fiTEx- 10 
(FEA) 

43.12 

18.15 

5.59 

3.92 

0.242 

0.222 

7.40 

25.44 

rule of 
mixt./ 
Halpin- 
Tsai eqns. 

43.40 

14.79 

4.45 

5.91 

0.260 

0.252 

6.77 

34.24 
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Table 17. 


Continuum 


properties for examples 5, 6 and 7 using SAM. 




E^Ct Ey 

(GPa) 

Ez 

(GPa) 

Gyz 

(GPa) 

Gxy 

(GPa) 

Vxz, 

v yz 

Vxy 

a x c ,a y c x 

lO^C 

a z c x 

10^/°C 

Example 5 

(plain- 

weave) 

uTEx- 20 
(SAM) 

12.46 

6.62 

1.64 

1.67 

0.399 

0.162 

29.10 

68.48 

fiTEx-10 

(FEA) 

11.81 

6.14 

1.84 

2.15 

0.408 

0.181 

28.36 

79.57 

Dasgupta 

results 

14.38 

6.25 

1.94 

3.94 

0.463 

0.167 

22.50 

86.00 

Example 6 

(plain- 

weave) 

fiTEx- 20 
(SAM) 

63.41 

11.13 

3.79 

4.24 

0.402 

0.027 

1.36 

21.53 

fiTEx-10 

(FEA) 

53.61 

10.88 

4.41 

4.72 

0.365 

0.128 

1.56 

22.7 1 

TEXCAD 

64.38 

11.49 

5.64 

4.87 

0.396 

0.027 

1.33 

20.71 

Example 7 

(5-hamess 

weave) 

fiTEx-20 

(SAM) 

69.30 

11.62 

4.06 

4.73 

0.355 

0.031 

1.21 

20.25 

fiTEx- 10 
(FEA) 

64.51 

11.33 

4.45 

4.85 

0.329 

0.047 

1.55 

22.03 

TEXCAD 

66.33 

11.51 

4.93 

4.89 

0.342 

0.034 

1.46 

21.24 


3.3 Selective Averaging Method (SAM) for Plate Properties 

In this section, the SAM procedure to compute the plate stiffness coefficients and 
plate thermal coefficients (for a thin textile composite) is described. To distinguish between 
the macrolevel, mesolevel and microlevel [A], [5] and fD] matrices, an over-tilde is used to 
denote the mesolevel stiffness, and a superscript ”M” is used to denote the macrolevel plate 
stiffness. However in the remaining sections, [A], [5] and [Z>] (without a superscript or 
overscript) will refer to the macroscale plate stiffness matrices. Also in this section, the 
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complete plate stiffness matrix on the macroscale, as defined by Eqn. (23) will be denoted 
by [C M l such that: 


[ C M ] 


A m B m 

B m D m 


(44) 


The procedure to find the plate stiffness and thermal coefficients is analogous to that 
used to find the continuum thermo-elastic constants. To find the first column of the effective 
plate stiffness matrix, the unit-cell is discretized into slices (mesolevel) and elements 
(microlevel) as shown in Fig. 20. The unit-cell is subject to the deformation given by eJJ=l. 
The following assumptions are made regarding the midplane strains and curvatures on the 
macrolevel, mesolevel and microlevel: 


= «.o = «« = 0 1 = 2,3 

xf = Xj = — 0 i — 1,3 


(45) 


It is also assumed that the non-zero strain component e x o, and the force resultant N x are 
uniform within the mesoscale and macroscale respectively. These two assumptions can be 
expressed as the following equations: 


£ *0 “ £jcO 

N x = Af 


(46) 


The mesolevel stiffness coefficient A { j can then be obtained by averaging the corresponding 
element stiffness coefficients over the slice (consequence of the isostrain assumption) as: 

c b 

A n (x) = \ | j Qn(x,y t z) dy dz (47) 

z=0y =0 

where Q \ i is the plane-stress stiffness coefficient in the classical lamination theory (Agarwal 
and Broutman, 1990), which has been transformed to the unit-cell’s xyz-coordinates (for an 
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isotropic material, Q u = 


E 


). The macroscale force and moment resultants can be 


11 (1 - V 2 ) 

expressed in terms of the microscale stresses by the following relations: 

c b a 

dx dy dz 


N; 


M 


~ ah j J J ff ' 


z = 0^ = 0jc = 0 
c b a 


i = 1,2,3 


(48) 


* r r 

z o t 

J J J 


~ \ I | z o i dx dy dz 


z = 0;y = 0;c=0 


The assumption of uniform force resultants on the macroscale and Eqn. (48) yields the 
following expressions for the first column of plate stiffness coefficients: 


a 

J_ = I f_J_ A 

A n J A n (x) 

x—0 


(49) 


c b a 


C n=ab j j j ^ Gi|(w) * dy dz (i 


= 1,2,3) 


(50) 


z=0y=0x=0 


f'M 1 

A ab 


c b a 


lih 

z=0y=0x=0 


A u (x) 


Qn(x,y,z ) dx dy dz 


(i = 1,2,3) 
j = i + 3 


(51) 


A similar procedure is followed to compute the second column of the effective plate stiffness 
matrix. 

In the case of shear loading = 1 , the unit-cell is discretized into slices parallel 
to the yz-plane. The assumptions of isostrain and uniform force resultants are reversed from 
the case of normal loading. The force resultant A^y is assumed to be uniform within a slice, 
such that 


Nxy = 


( 52 ) 
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It is also assumed that is the only non-zero deformation component on the macrolevel, 
mesolevel and microlevel. Averaging the element compliance coefficients over the slice, we 
obtain: 


b a 


— “ — = -7 I 7; — T 7 dx dy 

Q 33 {z) Ub J J ^33 (W) 


(53) 


_}’ = 0 x - 0 

The mesolevel stiffnesses are then averaged over the volume of the unit-cell to yield the third 
column of the plate stiffness matrix, as follows 

c b a 


r'M 

U 2 , 3 


• r r 

M _ 1 

i3 ah 

UD J J J 


Q i3 (x,y,z) 

Q 33 (x,y,z) 


Q 33 {z) dx dy dz ( i = 1,2,3) 


(54) 


z = 0;y = 0.x=0 


S3 


ab 


c b a 

111 

z=0y=0x=0 


Q i3 {x,y,z) 
Q 33 (x,y,z ) 


Q 33 (z) dx dy dz 


(i 

j 


1,2,3) 
i 4- 3 


(55) 


The expressions for the fourth, fifth and sixth columns of the plate stiffness matrix 
can be obtained in a similar fashion to that explained above, except that the one of the 
curvature component will be non-zero instead of a midplane strain component. For example, 
to determine the fourth column, x x will be the only non-zero deformation on the macroscale, 
mesoscale and microscale. Assuming that the curvature is uniform within a slice, we get 
x x ~ x x . Then by averaging the element stiffness coefficients over the slice, we obtain an 
expression analogous to Eqn. (47) as: 

c b 

f>ii(x) = - J J z 2 Q u (x>y,z) dy dz (56) 

z=0y=0 

The moment resultant M x is assumed to be uniform on the mesoscale such that M x = M x . 
By averaging the slice compliance coefficients over the unit-cell volume, we get the 
following relations for the fourth column of the plate stiffness matrix: 
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a 



( 57 ) 


c b a 


f f f D M 

C * ^ ab \ Z - 11 Qn( x >y*d d y dz (* = 

J J J D^(x) 


1,2,3) (58) 


z =0y=0x=0 


c b a 


r M 

c y4 


f [ f D M 

= ^J J J tj^Q^y’**** 


(i 

j 


1,2,3) 
/ + 3 


(59) 


z = 0y=0jc— 0 


The fifth and sixth columns of the plate stiffness matrix can be found by using a similar 
procedure. 

To find the plate thermal coefficients, c# and x t are assumed as zero in the 
macrolevel, mesolevel and microlevel. The thermal stresses developed in the microscale due 
to a uniform temperature difference of 417=70 are given by: 


{<>} = ~ [Q]{a}T 0 

where 



'Ox 


r a/ 

{<7} - - 

Oy 

^xy 

t - 

{a} - - 

a y 

0-xy 
•» - 


The macroscale plate constitutive equation will reduce to 


(60) 


(61) 



(62) 


By averaging the microscale stresses given by Eqn. (60) over the unit-cell volume (using 
Eqn. 48), and equating to the macroscale force and moment resultants in Eqn. (62), we get 
the following relations for the plate CTE’s: 




( 63 ) 
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3.4 Plate Results using SAM 


The plate coefficients and plate CTE’s were computed for examples 2-7 by 
implementing SAM, and compared with the finite element results presented in the previous 
chapter. The computed coefficients for the examples are listed in Tables. 18-20. The 
properties for the bimaterial plate were predicted very accurately (Table 18). For the 
unidirectional composite examples, SAM was able to predict the [A] matrix coefficients 
(except for A^) t D \\ j and the plate CTE’s very well (Table 19). However SAM grossly 
over-predicted D \ 2 and D 22 and under-predicted the coefficient D For the textile 
composite examples (Table 20) SAM predicted all but A 12 and D \2 with very good accuracy. 
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Table 18. Non-zero [A], [5] and [D] coefficients for example 2 (bimaterial plate) using 
SAM. 



An, A22 
xlO 6 

A ]2 

xlO 6 

A66 

xlO 6 

Blh B22 
xlO 3 

Bn 

xlO 3 

B 66 

XlO” 3 

uTEy- 20 
(SAM) 

9.844 

2.045 

3.899 

-0.565 

-0.108 

- 0.228 

fiTEx- 10 
(FEA) 

9.832 

2.043 

3.895 

- 0.563 

-0.108 

-0.228 

lamination theory 
for two plies 

9.832 

2.043 

3.895 

- 0.563 

-0.108 

- 0.228 

lamination theory 
using 3-D elastic 
constants 

9.844 

2.048 

3.899 

0 

0 

0 



Du, D 22 
xlO -3 

D 12 

xlO -3 

A % 
xlO -3 

a x P, a y P 
xlO- 6 /°C 

MV 

PC / m 

juTEy -20 
(SAM) 

53.727 

11.163 

21.282 

17.828 

0.170 

fiTEx - 10 
(FEA) 

53.590 

11.149 

21.220 

17.800 

0.170 

lamination theory 
for two plies 

53.573 

11.131 

21.220 

17.814 

0.170 

lamination theory 
using 3-D elastic 
constants 

53.762 

11.183 

21.293 

8.190 

0 


Note: [A], [5] and [D] coefficients in SI units 
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Table 19. Non-zero [A], [5] and [£>] coefficients for examples 3 and 4 (unidirectional 
composite) using SAM. 


An 

xlO 6 


Example 3 

uTEy- 20 
(SAM) 

0.688 

0.175 

0.475 

0.102 


liTE%- 10 
(FEA) 

0.690 

0.149 

0.496 

0.177 


Halpin-Tsai Eqns. 
and lamination 
theory 

0.673 

0.109 

0.363 

0.113 

Example 4 

uTEy- 20 
(SAM) 

0.443 

0.074 

0.261 

0.040 


ftTEy-M 

(FEA) 

0.452 

0.062 

0.285 

0.114 


Halpin-Tsai Eqns. 
and lamination 
theory 

0.444 

0.039 

0.151 

0.045 



Dn D 12 D 22 

xlO -6 xlO -6 xlO -6 



a x Px a y Px 

10-6 10 -6 


Example 3 

uTEy- 20 
(SAM) 

3.750 

1.029 

3.112 

0.541 

14.476 

24.750 


luTEx- 10 
(FEA) 

3.589 

0.596 

1.980 

0.947 

15.489 

26.184 


Halpin-Tsai 
Eqns. and 
lamination 
theory 

5.606 

0.908 

3.026 

0.939 

15.625 

55.112 

Example 4 

uTEx~ 20 
(SAM) 

2.308 

0.446 

1.799 

0.195 

6.628 

13.076 


liTEx- 10 
(FEA) 

2.256 

0.224 

0.873 

0.568 

7.378 

13.188 


Halpin-Tsai 
Eqns. and 
lamination 
theory 

3.702 

0.328 

1.262 

0.371 

6.774 

34.239 


Note: [A], [5] and [D] coefficients in SI units 
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Table 20. Non-zero [A], [5] and [£>] coefficients for examples 5, 6 and 7 using SAM. 


Example 5 

fiTEx-20 (SAM) 

2.667 

0.446 

0.379 

0 


fiTE%- 10 (FEA) 

2.681 

0.565 

0.489 

0 


lamination theory 
using 3-D constants 

2.783 

0.503 

0.490 

0 

Example 6 

fiTEx-20 (SAM) 

12.215 

0.577 

1.095 

0 


juTEx-10 (FEA) 

12.090 

3.470 

1.223 

0 


lamination theory 
using 3-D constants 

13.938 

1.787 

1.208 

0 

Example 7 

luTEx-20 (SAM) 

16.039 

0.631 

1.209 

0.590* 


fxTEx- 10 (FEA) 

14.683 

1.351 

1.210 

0.495* 


lamination theory 
using 3-D constants 

16.531 

0.770 

1.239 

0 


D]]> D 22 

A 12 

B>66 

a x P, CLyP 

A" 

x 10“ 3 

x 10“ 3 

x 10" 3 

x 10 -Vc 

/°C/m 


Example 5 

fxTEx-20 (SAM) 

6.017 

1.590 

1.360 

27.505 

0 


[iTEx- 10 (FEA) 

5.687 

1.518 

1.577 

27.465 

0 


lamination theory 
using 3-D 
constants 

12.054 

2.177 

2.124 

28.363 

0 

Example 6 

fiTEx- 20 (SAM) 

44.782 

2.733 

4.398 

1.984 

0 


fiTEx-10 (FEA) 

41.695 

0.373 

5.879 

1.480 

0 


lamination theory 
using 3-D 
constants 

75.942 

9.734 

6.582 

1.556 

0 

Example 7 

UTEx- 20 (SAM) 

81.164 

3.307 

5.868 

2.422 

-0.028* 


fiTEx- 10 (FEA) 

90.072 

1.123 

6.149 

2.910 

-0.037* 


lamination theory 
using 3-D 
constants 

87.283 

4.195 

6.753 

1.550 

0 


* In example 7, B 22 - -Bjj and pyP = -fi x P 

Note: [A], [5] and [D] coefficients in 


SI units 


































































































CHAPTER 4 

FINITE ELEMENT MODELS FOR STRENGTH PROPERTIES 


In Chapter 2, we had demonstrated finite element procedures to model a general 
textile composite either as a three-dimensional (continuum) material or as a thin plate^eam 
to predict their corresponding thermo-mechanical coefficients. In this chapter, we extend the 
same numerical models to compute the thermal residual stresses due to processing in the 
yams and the matrix. Then the numerical models are used to study the strength behavior of 
the composite by predicting the failure envelopes for thin and thick textile composites. 

4. 1 Thermally Induced Residual Microstresses 

The thermal residual microstresses are induced in the yam and matrix materials due 
to the mismatch in their corresponding CTE’s. The difference between the composite curing 
temperature and room temperature then serves as the driving force to create these 
microstresses. Since composites designed for high temperature applications are fabricated 
at higher temperatures, the residual microstresses become relevant in the strength 
considerations of such composites. The residual microstresses in the vicinity of the 
yam-matrix interface are particularly important as they could lead to failure due to 
debonding. 

Determination of residual microstresses . Let T 0 be the difference between room 
temperature and the composite fabrication temperature. Since the composite is stress free at 
the fabrication temperature, which is above room temperature, T a is generally negative. The 
residual microstresses in the yam and the matrix are obtained by superposing the 
microstresses due to the two load cases as explained below. In the first load case, the unit-cell 
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is constrained from expanding by fixing the comer nodes of the unit-cell and enforcing zero 
displacement difference between corresponding nodes on opposite faces of the unit-cell 
(periodic displacement BC’s). A temperature difference T 0 is applied to all elements in the 
finite element model. This is exactly the same problem we solve for finding the 
three-dimensional (continuum) CTE’s. The applied boundary conditions mean that all the 
macrostrain components are equal to zero ( {e M ] = 0 ,A T M = T 0 ). Then the corresponding 
macrostresses required to restrain the unit-cell expansion are given by: 

(ct m ) = - [C] {a c )T„ (65) 

In the second load case, deformations are applied so as to reverse the macrostresses 
developed in the first load case. This can be accomplished by imposing the deformations 
{e M } = { a c } T 0 and AT M =0. It can be noted that the macrostresses developed in the second 
loading case, given by [C\{a c )T 0 are equal and opposite to the macrostresses in Eqn. (65). 
The microstresses from both load cases are superposed to obtain the residual stresses due to 
free thermal expansion. 

The same idea can be extended to finding the residual microstresses in the plate 
model. Then the deformations to be applied in the first load case are { = 0, [x M ] - 0 and 
A T m = T 0 \ and the deformations in the second load case are { e q } = { a p) T q , {x M } = {j3 p } 
T a and A T M - 0. The residual microstresses were computed for the plain-weave textile beam 
at the Gaussian center of the elements in the unit-cell. The beam was assumed to be in the 
x-z plane with unit-cells repeating in the jc-direction. Figure 21 shows the thermal stress 
contours for o**, o zz and x xz . The composite curing temperature was assumed to be 150°C 
above room temperature. 
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4.2 Strength Modeling of Textile Composites 

There are many failure criteria or strength theories for unidirectional fiber 
composites. This, for example, includes maximum stress theory, maximum strain theory and 
Tsai-Hill theory (Agarwal and Broutman, 1990). Even though failure of a material is a very 
complex phenomenon, engineering strength theories such as mentioned above have been 
found to be useful in design. The interpretation of strength values obtained from such 
theories are different for different materials. For example in metal matrix composites the 
failure envelope obtained using the above theories will correspond to the initial yield surface 
(Dvorak et al., 1973). In graphite/epoxy composites the failure theories can be used to predict 
fiber or matrix failure. In the present study our intent is to explore the possibility of 
developing such failure criteria for textile composites. 

4.2.1 Determination of Composite Failure Envelope 

Our approach is similar to that used by Dvorak et al. (1973). A state of homogeneous 
deformation, corresponding to each of the six macrostrain components, are independently 
applied to the unit-cell by imposing the boundary conditions explained in Section 2.1. For 
each case, the various stress components are computed in the elements in the 
unit-cell — typically at the element Gaussian integration points. These stresses will be 
referred to as microstresses. Assuming linear elastic behavior, the microstresses can be 
computed for any arbitrary combination of the macrostrain components. Since we know the 
macroscale elasticity matrix, we can find a relation between microstresses at a point and any 
arbitrary state of macrostress as: 

M = [F] [a M ] (66) 

[FI can be considered as a matrix of influence coefficients, which is evaluated at the 
integration points of all the elements in the unit-cell. We also assume that the failure behavior 
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of the matrix material and the yam is known. For instance, let the failure criterion of the 
matrix be of the form [ H] {cr} matrix = 1. Then the failure criterion for the composite is 
obtained from Eqn. (66) as [H] [F] { o M } = 1. The same idea also applies for the yarn. The 
textile composite is assumed to have failed if there is failure on the microscale in any one 
of the constituent materials — either matrix or the yam. By varying the macrostresses using 
a numerical simulation, failure envelopes can be obtained for the idealized homogeneous 
material (Fig. 22a). It might be noted that Eqn. (66) can be modified to include the thermal 
residual stress field in the unit-cell as 

Iff) = [F] ( 0 "} - |ff T )r o (67) 

where [o T ] is the matrix of thermal microstresses computed at the element integration point 
for a unit temperature difference. 



Figure 22. Composite failure envelopes. 

(a) continuum failure envelope in the space of macrostresses; 

(b) plate failure envelope in the space of force and moment resultants. 



69 


4.2.2 Effect of Stress Gradients on Strength 

The strength analysis for a three-dimensional composite can be extended for thin 
composites using the plate model. As mentioned in Section 2.2, the macrostresses will not 
be homogeneous thorough the thickness in such composites. Then the composite failure will 
be determined by the stress gradients through the thickness, represented by the averaged 
force resultants (N) and moment resultants (M). The composite failure criterion for the plate 
model will be of the form: 

ramljl = i (68) 

Thus the failure envelope of the composite will be in the six-dimensional space of the force 
resultants and the moment resultants (Fig. 22b). The above procedure was demonstrated 
using the beam model, for the case of a plain-weave textile composite. 

Failure envelope for textile composite beam . For a textile composite beam, the 
failure envelope is constructed in the space of the three force resultants P, M and V. The textile 
composite beam is assumed to be in the xz-plane with unit-cells repeating in the x-direction. 
The unit-cell is discretized with eight-node isoparametric plane strain finite elements. Three 
linearly independent deformations, as explained in Section 2.3, are applied to the unit-cell: 
(a) unit axial strain; (b) unit curvature such that transverse shear strain is zero; (c) unit 
transverse shear strain. For each deformation, from the finite element results, the 
microstresses a xx , a zz and x xz are computed at the Gaussian center of each element. The 
microstresses in the / th element for a combination of loads (deformations) are given by: 


'< 4 ' 



feAO 

fc 0 

<4 

> 

= [A 1 ] < 


r 

\ 



f 

^ ... 


( 69 ) 
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where the each column of [A*] corresponds to the microstresses for unit deformation. For 
example, the first column of [A^] represents the three microstress components in the z' th 

element for unit £q. We also know that 


fc 0 


P' 

< x M 

> = [S] i 

M > 

,.M 


V 





(70) 


where [S] is the beam compliance matrix. From Eqns. (69) and (70) we arrive at an 
expression analogous to Eqn. (66): 


" * ■» 
O'xx 

. 

'P' 

* Ozz 

T'XZ 
V* J 

> = [F] < 

M 

V 


(71) 


where [ F] — [A 1 ] [S]. The composite is assumed to have failed if there is failure in any one 
of the finite elements in the unit-cell. Matrix failure is determined by either maximum 
principal stress or von Mises criteria. Since the yam is assumed to be transversely isotropic, 
yam failure is determined by either maximum strain theory for a unidirectional composite 
or the Tsai-Wu criteria. 


Beam Failure Envelope Results . The strength properties used for the constituent 
materials in the beam are as follows: 

Yam: cr L r =1725 MPa, a L c =1366 MPa, o T r = 42 MPa, o T c =230 MPa, t lt =95 MPa 

Matrix: o 7=70 MPa, o c - 100 MPa 

where the superscripts ’ T and ’ C refer to the tensile and compressive strengths respectively. 
The failure envelopes were developed using two different sets of failure criteria. In the first 
case maximum principal stress criterion was used to determine the matrix failure and the 
maximum strain theory for unidirectional fiber composite was used for the yam. In the 
second case von Mises criterion was used for the matrix and the Tsai-Wu criteria for the yam. 
Both structural (beam) and continuum failure envelopes were developed. As was explained 
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earlier, in the beam model, periodic BC’s were imposed between corresponding nodes on 
the left and right ends of the unit-cell. In the continuum model, the periodic BC’s were in 
addition imposed between the top and bottom surfaces. The failure envelopes in Figures (23) 
through (26) were obtained using the first set of failure criteria. Figures (23) and (24) depict 
the structural failure envelope in the P-M space based on yarn and matrix failure 
respectively. Figures (25) and (26) illustrate similar continuum failure envelopes in the space 
of the macroscale normal stresses in the x- and z-directions. As expected, the envelope 
reduced in size with increasing shear force resultant. If we assume that the beam is made of 
a homogeneous but orthotropic material with properties as predicted by the continuum 
model, then one can derive structural failure envelopes from the continuum failure envelopes 
using simple beam theories. The derived structural failure envelopes are compared with that 
obtained from direct micromechanical analyses in Figures (27) and (28). One can note that 
the continuum failure criteria are very conservative for the case of a thin beam. Figures (29) 
and (30) show similar comparisons for the second set of results obtained using the quadratic 
failure criteria for the constituent materials. 



bending moment' (M) 


72 



axial force (P) X 10 3 


Figure 23. Beam failure envelope based on yam failure, (maximum strain theory) 



axial force (P) x jqj 

Figure 24. Beam failure envelope based on matrix failure, (maximum principal 
stress criteria) 
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xlO* 



Figure 25. Continuum failure envelope based on yam failure, (maximum strain 
theory) 


xlO 7 



Figure 26. Continuum failure envelope based on matrix failure, (maximum principal 
stress criteria) 
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Figure 30. Comparison of failure en' 
criteria) 






CHAPTER 5 

ISSUES IN MICROMECHANICAL MODELING 


The complex yam architectures in textile composites make traditional finite element 
modeling very difficult. Traditional finite element models assume that the material 
properties are constant or vary smoothly within an element. In the context of textile 
composites, it means that the yam and matrix materials are modeled by separate elements, 
with common nodes at the yam-matrix interface. This is indeed preferable because the 
stresses at the interface can be computed accurately. However meshing the interstitial matrix 
region becomes time-consuming and difficult. This region which is essentially a collection 
of multiply-connected matrix pockets, requires a very fine mesh to capture the pocket 
geometry. Dasgupta et al. (1990) and Whitcomb (1991) presented detailed finite element 
analyses of the unit-cell for plain-weave architectures. In this chapter, alternative finite 
element modeling techniques (Marrey and Sankar, 1994) are discussed, which are valid for 
any yam architecture. The methods were tested with simple two-dimensional examples, and 
compared with the displacement fields using traditional finite elements. Some aspects of 
finite element mesh generation, in the context of textile composites, are also discussed. 


5.1 Finite Element Modeling of the Unit-Cell 


Traditional finite elements use homogeneous elements, i.e., elements that are 
comprised of only one material. Two methods to circumvent traditional finite element 
modeling are described in the following section. The first uses incompatible elements, which 
employs homogeneous elements for the yam and matrix, though there is a node mismatch 
at the yam-matrix interface. The second method utilizes inhomogeneous elements to model 
the unit-cell. 
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5.1.1 Incompatible Elements 

The following method is a modification of previously published work on 
independently modeling substructures with finite elements (Ransom et al., 1993). The 
matrix and each of the yams are meshed independently with finite elements. In general, the 
nodes on the surface of the yarns will not coincide with the matrix nodes on the yam-matrix 
interface. In fact, most textile composite unit-cells due to multidirectional reinforcements, 
contain multiply-connected matrix pockets. These matrix regions require an extremely fine 
mesh to capture their geometry, whereas the yams can be modeled with a coarse mesh. This 
method has therefore two advantages — the first being the relative ease of generating 
independent meshes for the yams and the matrix. Secondly, the effective degrees of freedom 
in the numerical model are reduced, since the yams are modeled with a coarse mesh. 

Figure 31(a) shows an example of a rectangular unit-cell with an ellipsoidal 
inclusion. The matrix mesh is shown in Fig. 31(b) and the inclusion mesh is shown in Fig. 
31(c). The matrix degrees of freedom on the inclusion-matrix interface are denoted as qb m \ 
the matrix degrees of freedom in the interior as qi m and the inclusion degrees of freedom are 
denoted as q y . The total strain energy in the unit-cell is given by 

U = \[ its <1 ? J 

+ - q y T Ry - qfR f - qfR% 

where K y is the finite element stiffness for the inclusion; and K m is the stiffness for the matrix 
material, which is divided into four sub-matrices. R y and R m are the external load vectors in 
the inclusion and matrix mesh respectively. The matrix degrees of freedom on the boundary 
are eliminated using the transformation, 

= Tq y (73) 

where Tis a transformation matrix. By substituting for qb m in Eqn. (72) and minimizing the 
strain energy with respect to q y and qi m , we get the following relation: 


’V m XT™’ 
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Figure 3 1 . Example to illustrate incompatible elements. 

(a) problem to be modeled; (b) matrix mesh; (c) inclusion mesh. 


Equation (74) can be solved for qy and qi m , and qb m can be recovered from Eqn. (73). The 
displacement at any point in the unit-cell is obtained by identifying the inclusion or matrix 
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element enclosing the point, and interpolating the nodal displacements of that element. By 
repeating this procedure for several representative points, the displacement field within the 
unit-cell is obtained. 

The method was first tested by predicting the displacement field for the unit-cell of 
a unidirectional composite subject to uniaxial tension. Due to symmetry of the unit-cell, only 
one-quarter of the unit-cell was modeled. The boundary conditions imposed required the 
edges of the unit-cell to remain straight after deformation (Fig. 32a). The initial mesh for 
the problem is shown in Fig. 32(b). The node mismatch at the inclusion-matrix interface can 
be observed in the figure. Eight-node isoparametric elements were used to model the fiber 
and matrix — 24 elements were used for the fiber and 98 elements for the matrix. The 
deformed incompatible element configuration is shown in Fig. 32(c). The results were 
compared with conventional finite element results by plotting the displacement fields for 
both cases on a background 25x25 mesh (Fig. 32d). It was found that the incompatible 
element displacement field was identical to that using traditional elements. Only one set of 
displacements are visible in Fig. 32(d), because the displacement fields exactly overlapped 
each other. 

The procedure to determine the background mesh displacement field is illustrated in 
Fig. 33. For a given point in the background mesh (Fig. 33a), the corresponding incompatible 
element containing the point is determined (Fig. 33b). The displacement for the mesh point 
is obtained as the interpolation of nodal displacements of the found incompatible element 
(Fig. 33c). A similar procedure is followed to compute the background mesh displacements 
using conventional elements. The method was also tested by modeling the transverse shear 
of a unidirectional composite (Fig. 34). In this case, one-half of the unit-cell was modeled 
and the displacements are shown for a background 20x20 mesh. In both examples, it was 
found that the displacement fields using incompatible elements and conventional elements 
were identical. A discontinuity in strains in the vicinity of the fiber-matrix interface was 
observed in both the examples. 
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(a) (b) 



Figure 33. Calculating displacement field for the background mesh. 

(a) background mesh showing node whose displacements are to be 
calculated; (b) incompatible element mesh; (c) matrix element 
containing the background mesh node. 
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Figure 34. Unidirectional composite subjected to transverse shear. 

(a) boundary conditions; (b) initial mesh using incompatible elements; 
(c) deformed mesh; (e) deformed background mesh. 
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5.1.2 Inhomogeneous Elements ( Averaged Gaussian Integration") 

By inhomogeneous elements, we mean finite elements which are comprised of more 
than one material. Inhomogeneous elements in micromechanical analyses were studied in 
detail by Foye (1993). The advantage of using inhomogeneous elements is that mesh 
generation is very simple. For instance a rectangular domain can be discretized into uniform 
rectangular or triangular elements, without taking into consideration the constituent material 
geometry. The stiffness matrix of inhomogeneous elements represent smeared properties of 
the constituent phases determined by the numerical integration scheme. Thus the solution 
will be approximate in the interface regions. Foye has developed a modified method to 
evaluate the stiffness matrix of inhomogeneous elements called Replacement Elements, 
which predict better stresses than conventional inhomogeneous elements. The following 
details the procedure to use inhomogeneous elements for micromechanical analysis. 

The unit-cell is divided into uniform rectangular hexahedral elements as shown in 
Fig. 35. In general these elements will be inhomogeneous. The stiffness matrix ( K e ) for an 
inhomogeneous element is formulated as: 

K e = B t CB dV e 

Jv * 

+1 +1 +1 

B t CB Ul d£drjdt (75) 

-1 -1 -1 

N N N 

= X2X W i W J W k B T C(i,V.& B\J\ 

/= 1 y= 1 jfc=] 

where V s is the domain of the element, B is the strain-displacement transformation matrix, 
C is the elasticity matrix, N is the number of Gauss points used for integration, W is the 
Gaussian integration weight factor and \J\ is the determinant of the Jacobian. The material 
property at each Gauss point is determined, and the corresponding elasticity matrix 
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is used to perform the volume integration. The element stiffness matrix will thus represent 
the averaged properties of the constituent materials in that element. 
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The approach was implemented by writing a code to identify the material property 
of an arbitrary Gauss point. The inclusion geometry was defined by a coarse finite element 
mesh. If the Gauss point fell within any of the inclusion elements — the material property of 
the point was identified as that of inclusion. Otherwise the material property of the point was 
identified as matrix. The algorithm to determine the element belonging to the Gauss point 
is explained in the following section. The results obtained by using the approach to the 
uniaxial tension of a unidirectional composite are shown in Fig. 36(a). The displacements 
were in agreement with those obtained using homogeneous elements (Fig. 36b). 

The problem with inhomogeneous elements is that they cannot represent the jump 
in strains that can occur at the yam-matrix interface. In fact there are three strain components 
that can be discontinuous at the interface, but the corresponding stresses must be continuous. 
Such a behavior cannot be represented by inhomogeneous elements which assume a 
continuous strain fields within the element. This problem can be resolved by decomposing 
the displacement field into two parts: a displacement field qj that produces a strain field 
continuous everywhere in the unit-cell, and the second one <72 that has a strain discontinuity 
at the yam-matrix interface. The field <72 can be assumed to be such that the displacements 
are identically equal to zero everywhere in the matrix and at the interface, and exist only in 
the interior of the inclusions. Thus one can use inhomogeneous elements for solving the first 
set of displacements. The second set of displacements exists only in the inclusions, and they 
can be solved by discretizing only the inclusion. However the issue is determining the 
decomposition q=qj+q2 • The condition for the decomposition is that the jump in interfacial 
stresses should be equal and opposite in the two problems, since the interfacial stresses are 
continuous in the given problem. 
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Point location algorithm . In the context of a textile composite problem, the material 
property code of a point is obtained as follows. The yam (inclusion) volume is discretized 
into eight-node hexahedral isoparametric finite elements as shown in Fig. 37. The finite 
elements are shaped as triangular prisms (a face of the hexahedron is collapsed to an edge). 
An iterative algorithm is used to determine whether the point (whose material property is 
to be determined) is contained within the volume of a given element. If not, the procedure 
is repeated for all the elements in the yam. If the point does not belong to any of the yams, 
it is designated as a matrix point. If the point belongs to a yarn element, the local yam 
direction is computed as the direction of the line along the length of the prismatic element. 



Figure 37. Finite element mesh for a yam. 


The algorithm to determine if a given point is contained within the volume of a finite 
element is explained below. The procedure is explained for the two-dimensional case, and 
can be easily extended to a three-dimensional problem. 
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Figure 38. Finite element in: 

(a) xy-coordinate system (b) natural coordinate system. 


Let the point, whose material property is to be determined, be P with coordinates (x,y). It 
is required to determine whether the point P is contained within the area of the eight-node 
isoparametric element as shown in Fig. 38(a). The element is mapped into its natural 
coordinate system (Fig. 38b). Similarly, point P is mapped to F*in the element natural 
coordinate system. The element in the (|,j/) natural coordinate system is a square of length 
two units with edges parallel to the ^ 77 -axes, such that the center of the square O* coincides 
with the natural coordinate origin. This means that we need to find whether the point P * is 
contained within the given square. This is possible if we know the natural coordinates of 
point P*. However there does not exist a transformation from the xy-coordinate system to 
the natural coordinate system. The reverse transformation exists as given below: 
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where denotes the shape function, and (x^y,-) represents the ^-coordinates for the i th node. 

Therefore we use an iterative procedure to determine the natural coordinates of P*. 

We substitute £ =0 and ^=0 in Eqn. (76), and determine the ^-coordinates of point 
O (corresponding to point O* in the natural coordinate system). The point O will not, in 
general, coincide with P. Then the error, i.e., the difference in the x- and y-coordinates 
between points P and O is calculated as Ax and Ay respectively. Also since the x- and 
y-coordinates are functions of the natural coordinates, Ax andzly may be expressed as: 


^ + 

Ay = + i Ar > 


(77) 


Equation (77) may be rewritten in matrix form as follows: 







(79) 


where [J(^,tj)] is the Jacobian at (£ 77 ). The Jacobian is computed using the following 
equation: 
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Knowing the Jacobian, A% and At} can be evaluated from Eqn. (79) as: 
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The parameters At; and Arj are the corrections to be applied to the coordinates of point O* 
to determine the coordinates of point P*. As a check, the natural coordinates of point P*is 
then transformed to its xy-coordinates using Eqn. (76), and compared with the given 
coordinates of point P. The procedure is repeated if the difference in coordinates is above 
a prescribed tolerance. 


5.1.3 Periodic Boundary Conditions 


The unit-cell is the smallest volume element in the composite which is representative 
of the yam architecture. The yam architecture in the composite is generally periodic, and the 
unit-cell contains one repeat of the yam pattern. Thus the composite structure can be formed 
by assembling the unit-cell in all three dimensions. When the composite is subject to any 
arbitrary load, the microstresses and displacement gradients will be continuous across the 
faces of the unit-cell. Continuity of microstresses requires that the tractions be equal and 
opposite at corresponding points on opposite faces of the unit-cell. Also the displacements 
between corresponding nodes on opposite faces of the unit-cell will differ only by a constant. 


Multipoint constraint elements . The periodic boundary conditions (traction and 
displacement boundary conditions) can be implemented using multipoint constraint 
elements, which are based on the principle of Lagrange multipliers. The finite element 
formulation with multipoint constraint elements will be of the form: 
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where {A} is the matrix of Lagrange multipliers (Cook et al., 1989). However this method 
requires large storage for the stiffness matrix, due to the degrees of freedom contributed by 
the constraint elements (A*). Also the resulting stiffness matrix will not be positive definite, 
and the matrix bandwidth will be very large. 
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Transformation equations . The following section explains an alternative method to 
enforce the periodic boundary conditions. The method is explained for the two-dimensional 
case, and can easily be extended for the three-dimensional problem. Consider a finite 
element mesh of a unit-cell (Fig. 39a), where the periodic boundary conditions are to be 
implemented between the left and right edges, and the top and bottom edges. Assume that 
each node has only one degree of freedom (d.o.f). The nodal degrees of freedom on the left 
and bottom edges (except the prescribed d.o.f’s) are designated as ’’master” degrees of 
freedom q m , and those on the right and top edges as ’’slave” degrees of freedom q s . The nodal 
degrees of freedom in the mesh interior are also classified as master degrees of freedom. 
There may arise a situation where the d.o.f for the comer node at x/=0, *2 = 0 is not 
prescribed. In that case, it is designated as a master node with two slave nodes, whose 
coordinates are given by (x;=0, X 2 = Lq) and ( xj=Lj,X 2 = 0). Thus a master d.o.f may have 
zero, one or two slave d.o.f’s depending on its coordinates. 

The slave d.o.f’s are transformed to master d.o.f’s by the transformation: 
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q s = Aq m + A q s (83) 

where A is a transformation matrix such that = d c(i y, and d denotes the Kronecker delta. 

”c” is an integer function such that for the i th slave d.o.f, c(i) is the corresponding master 
d.o.f. Equation (83) can be rewritten as. 


'y' A.jj(q m )j + (Aq s )i 
j 

= + ( g4 ) 
j 

~ (Qm) c (i) "f (Aq s )i 

{A q s )i 'is the displacement difference to be imposed between the i th slave d.o.f (q s )i and its 
corresponding master d.o.f, (q m ) c ^y The strain energy in the unit-cell is evaluated as, 
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where X is the matrix of Lagrange multipliers; R s and R m represent the external force vectors 
for the slave and master d.o.f ’s. By minimizing the strain energy with respect to the variables 
q m , q s and X we get: 
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Eliminating q s and X from Eqn. (86), we arrive at the expression : 

K*q m - R* 


where 


K* = [K mm + A t K ss A + A T K sm + K ms A ] 
R* = R m + A t R s - [K ms + A T K ss \Aq s 
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Equation (87) can be solved for q m \ and q s recovered from Eqn. (83). Thus the number of 
unknowns in the linear equation solver are reduced, and also the resulting stiffness matrix 
( K *) is positive-definite. 


5.2 Mesh Generation 

Currently approximate methods such as averaged Gaussian integration. Selective 
Averaging Method and Isostrain Method are popular because of their ease of computation 
and requirement of minimum computer storage and time. However the advancements in 
computer hardware and computational technology will make it possible to use a large 
number of degrees of freedom for micromechanical analysis. Then there will be a need for 
mesh generation techniques for creating homogeneous elements in a unit-cell. In the 
following we describe two methods for discretizing the unit-cell with homogeneous 
elements. 

5.2.1 Node Migration Method 

The unit-cell is meshed with tetrahedral solid elements in an arbitrary fashion. This 
will be called the primitive mesh. The elements are identified as homogeneous or 
inhomogeneous. Then the nodes of the inhomogeneous elements are moved to the interface 
by using an heuristic algorithm. In each inhomogeneous element, the node closest to the 
interface is allowed to migrate to the interface first. After each cycle, some elements will 
transform into homogeneous elements. Then the process is repeated until all the elements 
become homogeneous. The mesh thus generated is called the intermediate mesh. 

The intermediate mesh will have some elements distorted due to node migration. 
This distortion can be removed by subjecting the mesh to an annealing process, by which 
the distortions concentrated near the interfaces are distributed among other elements also. 
An example of this method is depicted in Fig. 40. The purpose here is to discretize a square 
unit-cell containing two circular inclusions. The figure shows the initial mesh containing 


0.-2 
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uniform square elements, the intermediate mesh of homogeneous elements and the annealed 
mesh. 



Figure 40. Example problem for node migration method. 

(a) primitive mesh; (b) intermediate mesh; (c) annealed mesh. 
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5.2,2 Modified Node Migration Method 

In this method an initial mesh is created as in the previous method. The elements are 
divided into three groups: homogeneous matrix elements, homogeneous yam elements and 
inhomogeneous elements which contain both matrix and yam. Then the inhomogeneous 
elements are redesignated as matrix elements, thus leaving only two kinds of elements. We 
will call this as the intermediate mesh. Thus the finite element representation of the yam will 
be smaller than the actual yam, i.e., the yam mesh will be fully contained within the actual 
yam. Then the nodes on the yam-matrix interface in the intermediate mesh are allowed to 
migrate to the nearest point on the actual yam-matrix interface. As before a 
three-dimensional finite element program with fictitious material properties is used in this 
step to obtain an annealed mesh. An example of the final mesh for a square unit-cell with 
two inclusions in the shape of a quarter-circle is shown in Fig. 41. 





CHAPTER 6 

CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 


The previous chapters presented some of the aspects involved in the modeling of 
textile structural composites. Micromechanical analyses were developed to predict the 
three-dimensional elastic constants and CTE’s for textile composite materials. The effect of 
stress gradients in thin textile composites was highlighted, and consequently, an independent 
micromechanical analysis was developed for such composites. The issues involved in 
modeling a thin textile composite were demonstrated by first modeling the composite as a 
beam and then as a plate. Two codes called /iTE %- 10 and fiTE %- 20 were developed in 
FORTRAN 77, which implemented the finite element procedure presented in Chapter 2, and 
the SAM procedure presented in Chapter 3 respectively. The predicted macroscale 
thermo-elastic constants compared very well with existing models for textile composite 
materials. 

The finite element analysis and SAM procedures assume that the unit-cell of the 
composite is shaped as a rectangular hexahedron. Both analyses can be easily modified for 
other unit-cell geometries. For example, in the finite element procedure, the relations for the 
periodic boundary conditions will change based on the unit-cell geometry. However the 
assumption, that the unit-cell is rectangular, limited the textile composite examples for 
which the above codes can be implemented. For instance, the unit-cell for a braided 
composite is hexahedral with the included angle between the edges equal to the braid angle. 
The braid angle, in general, will not be equal to a right angle. Thus the effective composite 
properties could be computed only for woven architectures, where the unit-cell is a 
rectangular hexahedron. 
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A considerable amount of time was invested in developing mesh generation 
algorithms and alternative finite element procedures, in the context of modeling textile 
composites. Some of these ideas are presented in Chapter 5. The mesh generation algorithms 
such as the node migration method were very effective for two-dimensional examples. 
However for the three-dimensional case, some of the elements near the interface collapsed 
to create degenerate and skewed elements. Further, extremely small tetrahedral elements 
were required to capture the interstitial matrix geometry. Computer storage and CPU time 
limitations made the use of such a finite element mesh impractical. Initially multipoint 
constraint elements were used to impose periodic boundary conditions in 10. To 

reduce computer storage, the transformation method (for periodic boundary conditions) and 
the skyline solver (Bathe, 1982) were incorporated in the finite element code. Further 
reduction in computer storage can be realized by using sparse matrix solvers instead of the 
skyline solver. 

The averaged Gaussian integration technique was used to compute the 
inhomogeneous element stiffness matrix infiTEx-10. The inhomogeneous element averages 
the properties of the constituent materials in the element, and assumes that the strain 
distribution is continuous within the element. Consequently, the inhomogeneous element 
cannot capture the discontinuity in strains across the yam-matrix interface. Thus the 
corresponding microstresses computed in the vicinity of the interface will be inaccurate. The 
accuracy of the microstresses may be improved by refining the (inhomogeneous element) 
mesh, but that would greatly increase the degrees of freedom in the numerical model. This 
emphasizes the need for effective finite element preprocessing codes for meshing the 
unit-cell of the textile composite with homogenous elements. It will also be useful to 
interface the finite element results with suitable graphics software. This will enable the user 
to get a feel for the problem by directly visualizing the unit-cell deformed configurations, 
and the microstress/microstrain distributions within the unit-cell. 
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Chapter 4 presented failure envelope results for a plain-weave textile beam. Implicit 
to the procedure for determining the failure envelope, was the assumption that microscopic 
failure in the unit-cell translated to macroscopic failure of the composite. The procedure is 
simplistic and conservative, but yet of use to a structural designer. More work remains to be 
done in obtaining results for failure envelopes using the strength models for textile 
composite continuum and for a composite plate, as explained in Chapter 4. An issue which 
this study, due to time limitations, does not address is parametric analyses to study the effect 
of changing constituent material properties and fiber volume fraction on composite 
properties. Both the finite element analysis and SAM can be easily extended for other 
macroscopic composite properties such as thermal conductivities, electro-magnetic 
properties and so on. 
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1. Introduction 


This manual describes the codes fiTE%—10 and juTEx~20 9 which are based on 
micromechanical analysis procedures to predict the overall (macroscopic) thermo-elastic properties 
for a textile composite. The analysis assumes that the textile composite material is comprised of 
repeating unit-cells or representative volume elements (RVE’s). In other words, the unit-cell can 
be considered as a building block, which by assembling in all three dimensions would form the 
parent textile composite material as shown in Fig. 1(a). On the microscale (comparable to the 
dimensions of the unit-cell), the unit-cell consists of the reinforcing yams, and the interstitial matrix 
material (Fig. lb). However, on the macroscale (comparable to the dimensions of the structural 
component), the composite is assumed to be homogeneous. The object of the micromechanical 
analysis is to determine the macroscopic stiffness coefficients and coefficients of thermal expansion 
(CTE’s) of the textile composite from the yam and matrix material properties, and the yam geometry 
within the unit-cell. 




Figure 1. (a) Textile composite continuum; (b) unit-cell; (c) textile composite plate 
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1 .1 Constitutive Models 

The micromechanical analysis assumes that the unit-cells span the composite in all the three 
dimensions. This assumption is true in the case of a thick textile composite. However, there are many 
structural applications in which thin composites are used. In such cases, there will be fewer 
unit-cells in the thickness direction (Fig. lc). Then the composite structure is modeled as a plate to 
determine its corresponding thermo-elastic coefficients. Consequently, we have two constitutive 
models - for thick (three-dimensional) textile composites and textile composite plates. 


1 .1 .1 3-D Textile Composite 

Consider a 3-D textile composite with unit-cells repeating in all three dimensions. Let the 
unit— cell be a rectangular parallelepiped with edges parallel to the coordinate axes xj, X 2 and xj. On 
the macroscale the composite is assumed to be homogeneous and orthotropic and the composite 
behavior is characterized by the following constitutive relation: 
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where {cr} and {£ } are the macroscale stresses and strains respectively;^ T is a uniform temperature 
difference throughout the unit-cell; {a s } and [C] are the macroscale CTE’s and orthotropic elasticity 
matrix to be determined. 


1.1.2 Textile Composite Plate 


The textile composite plate is assumed to be in the xy-plane with unit-cells repeating in the 
x and y directions. The plate is assumed to be homogeneous on the macroscale and given by the 
constitutive relation: 
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where £io>YiO and x t are the midplane axial strain, shear strain and curvature; a p and fi p . are the plate 
thermal expansion and thermal bending coefficients; N x and JWj are the axial force and bending 
moment resultants respectively in the homogeneous plate. [A], [£] and [D] are the plate extensional 
stiffness, coupling stiffness and bending stiffness matrices respectively. 

1.2 Brief Background of fiTEx-10 and juTEx-20 

As was mentioned earlier, fiTEx~10 and / iTEx—20 are based on micromechanical 
algorithms to predict the effective stiffness and thermal coefficients for a textile composite. The 
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algorithms are briefly explained in this section. The concepts behind [iTE%— 10 are explained in 
detail in references 1-6. Reference 7 summarizes references 1-6, and also explains /jlTE%—20's 
computational procedure. 

1 .2.1 y uTEx-10 Theory: 

fiTEx—10 , pronounced as microtech- 10, is a finite-element code to predict the textile 
composite properties. In fact, this code can be used for any composite in which the inclusions follow 
a repeated pattern - such as a unidirectional composite with fibers in a square array. Both the 3-D 
(continuum) properties as well as the plate properties can be computed. 


3 



Figure 2. (a) unit-cell discretized into uniform brick elements (b) eight-node brick element 


The unit-cell is assumed to be a rectangular hexahedron and divided into regular-sized 
eight-node brick elements as shown in Figure 2. The number of divisions along the 1-, 2- and 
3-axes are specified by the user. Periodic displacement and traction boundary conditions are 
imposed between corresponding nodes on opposite faces of the unit-cell. The periodic boundary 
conditions (BC’s) ensure that the displacements and stresses are continuous across the faces of the 
unit-cell. It should be remembered that the elements are in general, inhomogeneous, i.e., consist of 
more than one material. The element stiffness matrix is obtained as the averaged stiffness of the 
different materials in the element. The numerical procedure implemented in juTE%— 10 to compute 
the element stiffness matrix is explained in appendix A.l. 

To determine the stiffness coefficients, periodic displacement BC’s are applied such that only 
one of the macrostrain components in Eqn. (1) is nonzero. For example, to impose unit eh, the face 
xi=Li is pulled by a distance of L\ in the 1-direction relative to its opposite face given by xi=0. 
Simultaneously, the difference in displacements in the 2- and 3-directions for the above two faces; 
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and the displacement difference in all three directions for the remaining two pairs of faces are set 
to zero. From the finite element results, the macrostresses required to create such a deformation state 
are computed. Substituting the macrostresses in Eqn. (1) and knowing that en is the only nonzero 
strain component, the stiffness coefficients corresponding to the first column of [C] are obtained. 
A similar procedure is followed to compute the remaining five columns of [C], and the composite 
CTE’s. 


The microstresses for the six deformation cases, and the thermal residual microstresses for 
a given temperature difference (between composite curing temperature and room temperature) can 
also be predicted. However, the computed microstresses will be inaccurate in the vicinity of the 
yam-matrix interface. 

1.2.2 fiTEx~20 Theory: 

The code juTEx—20 implements the Selective Averaging Method (SAM) to estimate the 
composite thermo-mechanical coefficients. SAM is an approximate analytical method which is 
faster and easier to implement than juTEx~10. The unit-cell is discretized into slices on the 
mesoscale as shown in Fig. 3. The direction in which the unit— cell is sliced is dependent on the state 
of loading. For normal loads (first three load cases) the slices are perpendicular to the direction of 
loading. In the case of shear loading (fourth, fifth and sixth load cases), the slices are parallel to the 
plane of shear. The slices are in turn sub-divided into elements on the microscale. The slice stiffness 
and then the composite (macroscale) stiffness are determined by utilizing a combination of isostress 
and isostrain assumptions on the mesolevel and macrolevel respectively. In other words, the 
composite stiffness coefficients are obtained by selectively averaging the stiffness and compliance 
coefficients of the elements and the slices. It should be noted that the periodic BC’s are not enforced 
in SAM. For further explanation the user should refer to reference 7. 



Figure 3. Hierarchy of discretization for a unit-cell to implement SAM: (a) unit-cell; (b) slice; 
(c) element 
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2. Preparing Input Files 

This section of the manual describes the procedure to prepare input files for both pTE%—10 
and liTE%- 20. The input files are identical for both the codes, except that post-processing options 
are available in p,TE%—10. In both codes, the XYZ-coordinate system is used for the unit-cell, and 
the 123-coordinate system for the yams. The XYZ-coordinate system is also referred to as the 
global coordinate system (GCS), and the 123-coordinate system as the local coordinate system 
(LCS). 

2.1 Input for pTE%-10\ 

Essentially, there are two input files to pTE%—10 - textile.inp and yarn.inp. The file 
textile, inp reads the unit-cell geometry, the yam and matrix material properties, and the user options. 
Whereas file yarn.inp reads the yarn/s architecture within the unit-cell. The procedures to create the 
above two input files are described in the following sub-sections. The various steps to create the 
input are stated along with a short explanation. After all the steps are discussed, examples of 
complete input files are documented. 

2.1 .1 Steps to create textile.inp : 

The below stated steps (with the exception of the first step) consist of a character line, which 
may be followed by numeric data. The character line may be entered completely in upper-case or 
completely in lower-case. The numeric data is represented by integer/real variables. The text to be 
actually typed in the files are preceded and followed by dashed-lines, to differentiate them from the 
rest of the text. 

1. The first step is to enter the comment line/s. Each line can be up to 80 characters long, and 
multiple lines (up to 100) are allowed. A blank line should be used to denote the end of the comment 
lines, for eg.. 


textile.inp - plain weave textile composite plate example 
Thermal residual micro stresses are to be computed 


2. Enter the material properties of the yams (each yam can have different properties). This is 
done by assigning a number to every yam. Then the data is input as: 


yarn 

NYR 

1 , EL , ET , GLT , PNULT , PNUTS , ALPHAL , ALPHAT 
2 , EL, ET, GLT , PNULT , PNUTS , ALPHAL , ALPHAT 


NYR , EL, ET, GLT , PNULT , PNUTS , ALPHAL , ALPHAT 


NYR is an integer representing the number of yams. Since the yarns are assumed to be 
transversely isotropic, five elastic-stiffness properties and two thermal expansion coefficients are 
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required for each yarn. In the coordinate system shown in Fig. 4, the L-axis corresponds to the axis 
of the yarn. The T- and S-axes are mutually perpendicular axes on the transverse plane, which is 
also the plane of isotropy. EL and ET are the elastic moduli of the yam in the longitudinal and 
transverse directions respectively. GLT is the shear modulus in the L-T plane; PNULT is the major 
Poisson ratio in the L-T plane, PNUTS is the Poisson ratio in the T-S plane. ALPHAL and ALPHAT 
are the yarn CTE’s in the longitudinal direction and transverse directions respectively. 




L 


Figure 4. Cross-section of yarn showing the longitudinal (L) and transverse (T and S) 
directions. T-S plane is the plane of isotropy. 

3. Enter the matrix (isotropic) properties. 


matrix 

EM, PNUM, ALPHAM 


EM and PNUM are the elastic modulus and the Poisson ratio of the matrix material; ALPHAM is 
the coefficient of thermal expansion of the matrix. 

4. Define the bounds of the unit-cell as follows: 


uc-domain 
LX, LY, LZ 


where LX, LY and LZ are the lengths of the unit-cell in the three dimensions. 

5. Specify the unit-cell type, and the number of elements that the unit-cell is to be discretized. 
If the composite 3-D continuum properties are to be determined, then the input would be of the form: 


solid 
NX , NY , NZ 
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On the other hand, if the plate properties are to be computed , then the input would be of the form: 


plate 
NX , NY , NZ 


NX, NY and NZ are the number of divisions along the the three dimensions. Consequently, the total 
number of elements generated would be NXxNYxNZ. Both ’solid’ and ’plate’ options cannot be 
used simultaneously. 

6. This step is optional. If the following line is typed, 


forces 


then the average stresses in the case of ’solid’ option, or the average force and moment resultants 
in the case of ’plate’ option, acting on each face of the unit-cell are computed and printed in file 
textile. out This option should be used only if required, since it increases the code run-time. 

7. This step is also optional. If the microstresses within the unit-cell are to be computed then 
the following line is entered: 


stresses 


Then the microstresses are computed at the Gaussian center of each element, for the six deformation 
cases and printed out in six different files, namely, stressl.dat, stress2.dat, stress3.dat, stress4.dat, 
stress5.dat and stress6.dat (Appendix A.2). If the thermal residual microstresses are also to be 
computed, then the input would be 


tstresses 

DELTAT 


where DELTAT is the difference between the room temperature and composite fabrication 
temperature (in most cases DELTAT is negative). Then, in addition to the above six files, the file 
tsress.dat will be generated, containing the thermal microstresses. 

The commands to compute the microstrains are analogous to the above two commands, and 
are given by: 


strains 


or 


tstrains 

DELTAT 


Then the files which would contain the microstrain data are strainl.dat, strain2.dat, strain3.dat, 
strain4.dat, strain5.dat, strain6.dat and tstrain.dat (Appendix A.2). 
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Note: Only one of the above four commands should be entered in textile, inp. 
8. This is the final step, which indicates the end of the input file. 


finish 


The order of steps 5, 6 and 7 may be interchanged. However, the user is advised to follow 
the order of the listed steps. Steps 1-5 and Step 8 are sufficient for computing the composite 
thermo-mechanical coefficients. Steps 6 and 7 should only be used when additional post-processing 
information is required, i.e, average forces and microstress/microstrain data. 

2.1.2 textile. inp example 


Example of textile. inp for a plain-weave textile composite. 

yarn 

4 

1,144. 8d9 ,11. 73 d9 ,5. 52d9 ,0.23,0. 30,-0. 324d~6 , 14 . Od-6 
2,159. 8d9 , 10 . 90d9 , 6.40d9, 0.38, 0.30, 0.0450d-6,20.2d-6 
3 , 144 . 8d9 , 11 . 73d9 , 5 . 52d9 , 0.23,0.30,-0 . 324d-6 , 14 . Od-6 
4,159. 8d9 , 10 . 90d9 , 6. 40d9, 0.38, 0.30,0. 045 Od-6 , 2 0 . 2d-6 
matrix 

3 . 45d9 , 0 . 35d0 , 40 . Od-6 
uc-domain 

2 . 822d-3 , 2 . 822d-3 , 0 . 2557d-3 

solid 

10,10,7 

forces 

tstrains 

-150. OdO 

finish 


In the above example, the unit-cell contains four yarns. Yams 1 and 3 have identical material 
properties and so do yams 2 and 4. Since the ’solid’ option is used, the continuum properties for the 
composite are computed and printed into file textile.out. The average stresses on the faces of the 
unit-cell are also printed in file textile.out. The microstrains for the six deformation cases and the 
thermal microstrains (for a temperature difference of -150 units) are also determined. 

2.1.3 Steps to create yarn.inp : 

The yam geometry information is entered in file yarn.inp. It is assumed that the end-faces 
of the yam are planar which coincide with the faces of the unit-cell. The yams are modeled by 
defining the end cross-section, and by sweeping this cross-section along the yam axial direction (or 
yam path). Listed below are the steps required to create yarn.inp. 

1 . Enter the comment line/s. A blank line should follow the last comment line, for eg.. 



Ill 


yarn.inp - plain weave textile composite plate 


example 


2. Enter the number of yarns as shown below. 


HEADERLINE 

NYR 


The HEADERLINE can be any string of characters (up to 80 characters), or even a blank line. It 
is essentially a dummy line, and is not used in the code. NYR stands for the total number of yams, 
and this integer must be the same as the total number of yams entered in textile, inp. An example of 
this step would be: 


number of yarns 
10 


3. Enter the number of points representing the cross-section of yam 1 as: 


HEADERLINE 

NPPTS 


where NPPTS denotes the number of points. For example, eight points are chosen in Fig. 5(a) to 
model the yam cross-section. Hence the input may be entered as: 


number of polygon points in yarn 1 
8 


4. Enter the (X 3 , X 2 ) coordinates for all the points representing the cross-section of yam 1 . This 

is done by choosing any point as the starting point and numbering all the points (from 1 to NPPTS) 
along the counter-clockwise direction (denoted by circled numbers in Fig. 5a). Then the (X 3 , X 2 ) 
coordinates are entered for every point in ascending order, as shown below: 


HEADERLINE 
X3 1 , X21 
X32 , X22 
X33 , X23 


X3 NPPTS , X2NPPTS 
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The first digit after ’X’ refers to the coordinate axes number and the second to the point number. For 
example X32,X22 would be the (X3, X2) coordinates of point 2 on the cross-section. The yam 
cross-section will therefore be approximated as a NPPTS-sided polygon as shown in Figure 5(b). 

5. Enter the number of points chosen to represent the axial path of yam 1. In the illustrated 
example 11 points are chosen to model the yarn path (Fig. 5c). The information is typed as: 


HEADERLINE 

NAPTS 


where NAPTS is the total number of axial points. 



2 



(c) (d) 

2 



2 



Figure 5. (a) yam cross-section in the 2-3 plane showing the points chosen to represent the 
yam cross-section; (b) approximated polygonal cross-section; (c) yarn path in the 1-2 plane; 
(d) approximated yam path. 


6. Input the (xi , X 2 ) coordinates for the axial points as follows (it is assumed that the yam path 
lies on the 1-2 plane) . First, number the points successively along the yam direction (Fig. 5c). The 
starting point should coincide with the origin. Then enter the data as: 
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HEADERLINE 
XI 1 , X21 
X12 , X22 


X INAPT S , X2NAPTS 

The first digit after ’X’ refers to the coordinate axes number and the second to the axial point number. 
For example X12,X22 are the (xi, X 2 ) coordinates of axial point numbered 2. Like in Step 4 the data 
should be entered in ascending order such that the last line corresponds to the coordinates for point 
NAPTS. 

7. Enter the transformation matrix for yam 1 to transform the coordinates from the yam 123 
coordinate system (LCS) to the global XYZ coordinate system (GCS). The transformation is given 
by the following equation: 

XI ^11^12^13 ( x \ 

< Y > = T 21 T 2 2 T 2 i < x 2 > (3) 

zj T 2 1 T32 T33 *3 

where [T] is the 3x3 transformation matrix. 


Z 



X 


Figure 6. Sketch to illustrate the transformation of coordinates from the yam coordinate 
system to the global cooridnate system 
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Consider Fig. 6, in which the yam coordinate system has been translated so that its origin 
coinicides with the global coordinate system. The components of the transformation matrix are 
determined by using the following equation: 

'(Xxj) (X,x 2 ) (X, x 3 y 

[T] = ( Y ’ X 0 ( Y * x 2> ( Y ’ x 3~) (4) 

(Z,x x ) (Z,x 3 ) (Z,x 3 ) 

where (Xi, xj) denotes the cosine of the angle between the axis Xj in the GCS and the jth axis in the 
LCS. Figure 7 shows an example of a unit-cell with one yam. Only the end-faces of the yam are 
shown. The 1-2 plane which is assumed to contain the yam path (not shown) is parallel to the Y-Z 
plane. From Eqn. (4) the transformation matrix for the yarn in the example would be: 

[0 0 - f 

[T] = 0 1 0 

L 1 0 °. 

The transformation matrix is entered as: 

HEADERLINE 
Til / T12 / T13 
T2 1 , T22 , T23 
T31 , T32 , T33 


Y 



Figure 7. Example to explain the transformation of coordinates from yam coordinate system 
to global coordinate system 
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8. Enter the global coordinates (xg,yg,zg) of the first axial point of yam 1, i.e., the global 
coordinates corresponding to the origin of the 123 coordinate system. This is typed as: 


HEADERLINE 
XG, YG, ZG 


9. Repeat Steps 3-8 for the remaining yams in the order of the yam number. 

2.1.4 yarn.inp example 

The following example is for a unit-cell with two yams. The unit-cell is of unit length in 
all three dimensions. The path of the first yam is assumed to be straight and parallel to the z-axis 
(Fig. 7). Therefore only two axial points are required to represent the path of the first yam. 


yarn.inp: Example of a composite with a two yarns. The first yarn 
is a straight yarn along z-direction with circular cross-section. 

number of yarns 
2 

number of points to represent cross-section of yarn 1 
12 

X3,X2 of polygon points 
0.437,0 

0.37845, 0.2185 
0.2185,0.37845 
0,0.437 

-0.2185,0.37845 

-0.37845,0.2185 

-0.437,0 

-0.37845,-0.2185 

-0.2185,-0.37845 

0,-0.437 

0.2185,-0.37845 

0.37845,-0.2185 

number of axial points to represent yarn path of yarn 1 
2 

XI, X2 of axial points 
0,0 
1,0 

transformation matrix for yarn 1 
0,0,-l 
0 , 1,0 
1 , 0,0 

global coordinates of first axial point for yarn 1 
0 . 5 , 0 . 5 , 0 

number of points to represent cross-section of yarn 2 
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4 

X3,X2 of polygon points 

0.4,0 

0,0.4 

-0.4,0 

0 ,- 0. 4 

number of axial points to represent yarn path of yarn 2 

5 

XI, X2 of axial points 

0,0 

0.25,0.41 

0.50,0 

0.75,-0.41 

1,0 

transformation matrix for yarn 2 
1 , 0,0 
0 , 1,0 
0 , 0,1 

global coordinates of first axial point for yarn 2 
0,0.5, 0.5 


2.2 Input for pTEx~20: 

The input files to pLTE%—20 are sam.inp and yam.inp (Appendix A.3). File sam.inp is 
similar in structure to that of textile, inp. The user should follow Steps 1-5 and Step 8 in sub-section 
2.1.1 to create sam.inp. Steps 6 and 7, which are post-processing options are not available for 
pTE%—20. The steps to generate yam.inp are described in sub-section 2.2.3. 
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3. Compiling the Codes 

Both juTEx-10 andfiTEx-20 are coded in FORTRAN. Double precision is used for all the 
real variables. The user might require to compile the codes depending on the size of the problem 
he/she needs to run. For instance, consider the case where the dimension statements for the number 
of divisions of the unit-cell in juTEx—10 were set to, say, seven in all three directions. If the user 
wanted to run a case with more than seven elements in any one of the directions, he/she will have 
to recompile p.TEx—10. Or conversely, the user might want to reduce the dimensions for a smaller 
problem, so as to reduce the executable size. In the following two sub-sections, the procedures to 
compile the codes are explained. 


3.1 Compiling pTEx~10 

The file paramlO.var contains PARAMETER statements which assign integer values to 
different parameters (or symbolic names). The parameters are in turn used to dimension the arrays 
in /iTEx—10. So the values assigned to the parameters are to be changed, to change the dimension 
of the arrays in /iTEx—10. 

Actually, only seven parameters in paramlO.var need to be changed, and the remaining 
parameters are automatically modified. The section of paramlO.var with the seven parameter 
statements is given below (for an example): 


c FILE PARAMl 0 . VAR 


c The variables in the following SEVEN 


-changed to 

change 

the array sizes 

parameter 

{NXDIM 

= 15) 

parameter 

{NYDIM 

= 15) 

parameter 

(NZDIM 

= 10) 

parameter 

(NSTIFF 

= 8341113) 

parameter 

(KYARN 

= 10) 

parameter 

(KAXIAL 

= 100) 

parameter 

(KPOLY 

= 100) 


parameter statements are to be 
the source code microtechlO. 


c the following dimensions are self-generated- 


parameter (NPOSX = NXDIM + 1) 
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The parameters NXDIM, NYDIM and NZDIM denote the number of divisions used to 
discretize the unit-cell along the X-, Y- and Z-directions. Fo example, if 17 elements are to be used 
along the X-direction, then the value assigned to NXDIM must be changed to 17 (or greater than 
17). KYARN stands for the number of yams in the unit-cell; KAXIAL represents the maximum 
number of axial points used for a yam; and KPOLY for the maximum number of cross-section points 
used for a yam. For instance, if 6 axial points are used to model one yam and 8 axial points for 
another, then KAXIAL must be greater than or equal to 8. 

The parameter NSTIFF assigns a value to the dimension of the global stiffness array in the 
source code. The global stiffness is stored as a one-dimensional array, using the skyline technique. 
In this technique, all the stiffness elements below the skyline 8 are stored in the one-dimensional 
array. The question arises, as to what value is to be assigned to NSTIFF. This is found using a 
two-step process. First enter the right values for the remaining six parameter statements, and assign 
a small number, say one, to NSTIFF. Then create the input files, textile.inp and yam.inp, compile 
and run the code. The following message should show up on the terminal, namely, 


total storage required for stiffness = NRSTIFF 

stiffness matrix is under-dimensioned 

set NSTIFF in file paramlO.var to required dimension 


where NRSTIFF is the exact storage required for the stiffness array. This integer value should be 
assigned to NSTIFF in paramlO.var and the code should be recompiled. The stiffness storage will 
depend on the number of elements used and whether the ’plate’ or ’solid’ options were chosen. The 
values to be assigned to the stiffness parameter for some discretizations (NXxNYxNZ) of the 
unit-cell are listed below. 


NX 

NY 

NZ 

NSTIFF (solid) 

NSTIFF (plate) 

5 

5 

5 

49173 

39948 

6 

6 

6 

126201 

94773 

7 

7 

7 

278001 

197298 

8 

8 

8 

549051 

373371 

9 

9 

9 

998733 

656832 

10 

10 

10 

1703493 

1090593 

11 

11 

11 

2759001 

1727718 

12 

12 

12 

4282311 

2632503 

13 

13 

13 

6414021 

3881556 

14 

14 

14 

9320433 

5564877 


For the above discretizations, the first step of the compilation procedure may be bypassed, i.e, the 
correct NSTIFF parameter can be assigned rightaway. An additional listing of stiffness parameters 
for discretizations ranging from 8-12 elements in the three directions are given in Appendix A.4. 

3.2 Compiling juTEx-20 

The parameters in file param20.var should be changed to alter the dimensions of arrays in 
code /uTEx—20. Given below is a partial listing of param20.var for an example. 
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c FILE PARAM2 0 . VAR 


c The variables in the following SIX parameter statements are to be 

c changed to change the array sizes in the source code microtech20. 


parameter (NXDIM = 50) 
parameter (NYDIM = 50) 
parameter (NZDIM = 50) 


parameter (KYARN = 10) 
parameter (KAXIAL = 100) 
parameter (KPOLY = 100) 


c the following dimensions are self-generated- 

parameter (NELEMS = NXDIM*NYDIM*NZDIM) 


Here NXDIM, NYDIM and NZDIM represent the dimensions for the number of divisions of the 
unit-cell in the X-Y- and Z-directions (Fig. 3). The definitions for KYARN, KAXIAL and 
KPOLY are the same as that in the previous sub-section. 
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4. Interpreting Ouput Files 

In this section, a brief explanation is given to understand the computed data for both the 
codes. Most of the data printed out in the output files is self-explanatory. Only the information which 
may not be obvious to the user is presented here. 

4.1 iiTE%-10 Output 

As shown in Appendix A.2, the output files are textile.out and yam.dat. In addition, if the 
’stress’ or ’strain’ options are used, then the stress/strain data are printed out in 
stress l.dat,...,stress6.dat or strainl.dat, strain6.dat (if the ’tstress’ or ’tstrain’ option were used 
instead, then an additional file tstress.dat or tstrain.dat is created). 

4.1.1 Interpreting textile.out 

Consider the following partial listing of file textile.out: 

plain-weave textile composite 
No. of Yarns = 4 


Yarn 

No . eL 

eT 

gLT 

pnu LT 

pnu TT 

1 

0 . 14480E+12 

0 . 1173 0E+11 

0 . 55200E+10 

0 . 23000E+00 

0 . 3 0000E+00 

2 

0 . 14480E+12 

0 .11730E+11 

0. 55200E+10 

0 . 23000E+00 

0 .30000E+00 

3 

0 . 14480E+12 

0 . 11730E+11 

0 . 55200E+10 

0 . 23000E+00 

0 .30000E+00 

4 

0 .14480E+12 

0 . 11730E+11 

0 . 55200E+10 

0 .23000E+00 

0 . 30000E+00 

Yarn 

No . alpha L 

alpha T 




1 

- . 32400E-06 

0 . 14000E-04 




2 

- . 32400E-06 

0 . 14000E-04 




3 

- . 32400E-06 

0 . 14000E-04 




4 

- . 32400E-06 

0 . 14000E-04 





Matrix Properties : 

E pnu alpha 

0 . 34500E+10 0 . 35000E+00 0.40000E-04 

x y z 


unit-cell length : 0.28220E-02 0.28220E-02 0.25570E-03 

divisions of SOLID unit-cell : 15 15 10 


volume 

fraction 

of 

matrix. . 

. 

= 0.14000 

volume 

fraction 

of 

yarn. . . 

1 

= 0.21506 

volume 

fraction 

of 

yarn. . . 

2 

= 0.21506 

volume 

fraction 

of 

yarn. . . 

3 

= 0.21500 

volume 

fraction 

of 

yarn. . . 

4 

= 0.21489 
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The volume fractions are computed by determining the material property code for all the element 
Gauss points. Eight Gauss points are used per element, which means that the total number of 
representative points will be eight times the number of elements. So the computed volume fraction 
indicates the fraction of the representative points whose material code corresponds to a particular 
material constituent. 

If the ’forces’ option were included in textile, inp, then the averaged stresses/force resultants 
are printed in textile. out. The nodal forces are computed at the end-faces of the unit-cell, and the 
nodal forces are averaged over the face to give the macroscale stresses or force resultants. For 
example, if the ’solid’ option is used along with ’forces’ option then the macroscale o xx is computed 
as: 

°xx = i Xr^ff(L x ,Y,Z) (5) 

^ n 

where F { f> is the nodal force in the x-direction at the nth node, and ^ denotes summation over all 

n 

nodes on the face X-Lx. On the other hand, for the ’plate’ option, the force resultant N x and moment 
resultant M x on the face X=Z* are computed as: 

n x = 

^ n 

1 v-. 

M x = f}^zF^{L x ,Y,Z) 

^ n 

Listed below are the macrostress data for an example problem (when ’solid’ and ’forces’ options 
were used): 


Results for SOLID MODEL obtained by AVERAGING FORCES 


Load case = 1 

Sigma XX 

Negative yz -plane 0.57691E+11 

Positive yz-plane 0.57691E+11 


Tau XY 

Negative xz-plane -.51854E+01 

Positive xz-plane -.51803E+01 


Tau XZ 

Negative xy-plane -.19312E+04 

Positive xy-plane -.19312E+04 


Tau XY 
- . 51828E+01 
- . 51828E+01 

Sigma YY 
0. 84979E+10 
0 . 84979E+10 

Tau YZ 
- . 59790E+04 
-. 59790E+04 


Tau XZ 
- . 19312E+04 
- . 19312E+04 

Tau YZ 
- . 59790E+04 
- . 59790E+04 

Sigma ZZ 

0 . 47023E+10 
0 . 47023E+10 


Load case 


2 
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3-D Continuum Stiffness Matrix is . . . 


0 . 57691E+11 

0 . 84979E+10 

0 . 84979E+10 

0.57674E+11 

0.47023E+10 

0.47020E+10 

- . 59790E+04 

- . 11463E+05 

- . 19312E+04 

- . 28005E+04 

- . 51828E+01 

- . 42448E+01 

Ex = 
Ey = 
Ez = 

0 . 55005E+11 
0 . 54988E+11 
0 . 10923E+11 

PNU yz = 
PNU xz - 
PNU xy = 

0 .35772E+00 
0 . 35774E+00 
0 . 11818E+00 

Gyz = 
Gxz = 
Gxy = 

0 . 43305E+10 
0.43308E+10 
0 . 4783 8E+10 


0 . 47023E+10 - . 59790E+04 
0 . 47 020E+10 - . 11463E+05 
0.11591E+11 - . 52929E+01 
-.52 93 0E+01 0 . 43305E+10 
- . 10377E+04 - . 46506E+01 
- . 93 648E+00 0.56016E+01 


- . 19312E+04 - . 51827E+01 
- . 28005E+04 - . 42448E+01 
- . 10377E+04 - . 93 648E+00 
- . 46505E+01 0 . 56016E+01 
0 . 43308E+10 - . 36642E+02 
- . 36642E+02 0.47838E+10 


3-D Coefficients of Thermal Expansion are . . . 

Alpha X = 0 . 15570E-05 
Alpha Y = 0 . 15580E-05 
Alpha Z = 0 . 22755E-04 

Alpha YZ = 0 . 15322E-10 
Alpha XZ = 0 . 80251E-11 
Alpha XY = 0 . 11073E-13 


The load case number refers to the state of deformation in the unit-cell. Load cases 1 to 6 refer to 
the state of deformation to determine the corresponding column of the plate/continuum stiffness 
matrix, and load case 7 refers to the state of thermal loading (to determine CTE’s). The data listed 
below the load case numbers are the averaged stress resultants. This means that the sign of the 
stresses on the negative planes (planes with outward normals in the negative X-,Y- and 
Z-directions) are reversed. The same applies to the force and moment resultants in case the ’plate’ 
option is used. These resultants obtained by averaging the nodal forces are used to construct the 
effective stiffness matrix - column by column. This stiffness matrix in theory has to be symmetric. 
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In practice, the finite element scheme and round-off errors may lead to unsymmetry of the stiffness 
matrix. Therefore, the effective compliance matrix is averaged, and the 3-D elastic constants are 
predicted by comparing the compliance coefficients with that of an orthotropic medium (this step 
is skipped for the ’plate’ model and the [A], [B], [D] coefficients are listed directly). 

The macroscopic stresses can also be obtained by averaging the corresponding microstress 
component over the unit-cell volume. For instance, the macroscopic cr** can be computed as: 

o xx = ±Jc%(X,Y,Z)dV (7) 

V 

where V is the volume of the unit-cell and er£ indicates the microstress component. In fact, 
irrespective of the options used , the effective stiffness matrix is always computed by averaging the 
microstresses (in addition, the stiffness is computed by averaging the nodal forces if the ’forces’ 
option is used). This stiffness matrix is listed at the bottom of textile.out as shown below: 

Results for SOLID MODEL obtained by AVERAGING MICROSTRESSES 
3-D Continuum Stiffness Matrix is ... 


0 . 57691E+11 

0 . 84979E+10 

0 . 47023E+10 

- . 59790E+04 

- . 19312E+04 

- . 51828E+01 

0 . 84979E+10 

0 . 57 674E+11 

0 . 47020E+10 

- . 11463E+05 

- . 28005E+04 

- . 42448E+01 

0 . 47 023E+10 

0 . 47020E+10 

0 . 11591E+11 

- . 52930E+01 

- . 10377E+04 

- . 93 648E+00 

- . 59790E+04 

- . 11463E+05 

- . 52930E+01 

0 . 43305E+10 

- . 46505E+01 

0 . 56016E+01 

- . 19312E+04 

- . 28005E+04 

- . 10377E+04 

- . 46505E+01 

0 . 43308E+10 

- . 36642E+02 

-. 51828E+01 

- . 42448E+01 

-. 93648E+00 

0 . 56016E+01 

- . 3 6642E+02 

0 .47838E+10 


Ex = 0 . 55005E+11 
Ey = 0 . 54988E+11 
Ez = 0.10923E+11 


PNU yz = 0 . 3 5772E+00 
PNU xz = 0.35774E+00 
PNU xy = 0 . 11818E+00 

Gyz = 0 . 43305E+10 
Gxz = 0 . 43308E+10 
Gxy = 0 . 47 83 8E+10 

3-D Coefficients of Thermal Expansion are . . . 



Alpha X = 0 . 15570E-05 
Alpha Y = 0 . 15580E-05 
Alpha Z = 0 . 22755E-04 
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Alpha YZ = 0 . 15322E-10 
Alpha XZ = 0 . 80251E-11 
Alpha XY = 0.11073E-13 


The method by which the elastic constants are calculated from the stiffness matrix is identical to the 
method explained in the previous paragraph. 

4.1 .2 Interpreting stress data 

If ’stress’, ’tstress’, ’strain’ or ’tstrain’ options are used, the microstresses will be printed out 
in different files as shown in Table 1. The first column states the deformation for which the 
microstresses are obtained. 


Table 1 . List of files containing the microstresses for various options. 



’stress’ option 

’tstress’ option 

’strain’ option 

’tstrain’ option 

1 

stress l.dat 

stressl.dat 

strainl.dat 

strainl.dat 

solid: Eyy=l 
plate: e yy =l 

stress2.dat 

stress2.dat 

strain2.dat 

strain2.dat 

solid: e 22 =l 
plate: y xy =l 

stress3.dat 

stress3.dat 

strain3.dat 

strain3.dat 

solid: Yyz=l 
plate: x xx =l 

stress4.dat 

stress4.dat 

strain4.dat 

strain4.dat 

solid: 7 x 2=1 
plate: Xyy=l 

stress5.dat 

stress5.dat 

strain5.dat 

strain5.dat 

solid: Yxy=l 
plate: x xy =l 

stress6.dat 

stress6.dat 

strain6.dat 

strain6.dat 

solid: thermal stresses 
plate: thermal stresses 

— 

tstress.dat 

— 

tstrain.dat 
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The following is a example listing of the file stress5.dat when the ’plate’ option was used: 

stress5 . dat 

micro stresses computed for deformation kappa-yy = 1 

The stresses are in the LTS-coordinate system for the yarns and 

the XYZ-coordinate system for the matrix. 

x y z yarn no. sigma LL sigma TT sigma SS 

tau TS tau LS tau LT 


0.47033E-03 0.47033E-03 
0.14110E-02 0 . 47033E-03 
0.23517E-02 0.47033E-03 
0.47033E-03 0.14110E-02 
0 . 14110E-02 0 . 14110E-02 


0 . 42617E-04 3 -.85259E+07 
0 . 18250E+04 
0 . 42617E-04 0 -.21669E+06 
- . 14835E+05 
0 . 42617E-04 1 -.11232E+06 
- . 20189E+06 
0 . 42617E-04 0 -.98276E+05 
0 . 71882E+04 
0 . 42617E-04 0 -.81038E+05 
0 . 42131E+03 


- . 16183E+06 0 . 90249E+05 
- . 50603E+05 - . 86194E+05 
- . 47759E+06 -.14582E+06 
0 . 10784E+04 0 . 33190E+05 
- . 16823E+07 - . 14407E+06 
0.99910E+04 0 . 57055E+05 
- . 24395E+06 -.33342E+05 
0 . 82 057E+03 0.32047E+05 
- . 21669E+06 - . 21813E+05 
- . 46568E+01 - . 64133E+05 


0. 23517 E- 02 0.23517E-02 0.21308E-03 2 0.11231E+06 0.16826E+07 0.14431E+06 

0 . 20124E+06 0 . 99994E+04 0.56985E+05 

Maximum tensile/compressive stress in each constituent material : 

yarn no. sigma LL sigma TT sigma SS tau TS tau LS tau LT 

0 0 . 21679E+06 0 . 47773E+06 0.14596E+06 0.14843E+05 0.10784E+04 0.33190E+05 
- . 21680E+06 - . 47773E+06 -.14597E+06 -.14835E+05 -.10958E+04 -.64133E+05 

1 0 . 15747E+06 0.16823E+07 0.14408E+06 0.13981E+06 0.10974E+05 0.27012E+06 
- . 15750E+06 - . 16823E+07 -.14407E+06 -.20189E+06 -.10984E+05 -.15633E+06 

2 0 . 15755E+06 0.16826E+07 0.14431E+06 0.20124E+06 0.10984E+05 0.59024E+05 
- . 15752E+06 - . 16826E+07 -.14432E+06 -.24088E+05 -.10993E+05 0.00000E+00 

3 0 . 85229E+07 0.16186E+06 0.90249E+05 0.47178E+05 0.49995E+05 0.13948E+06 
- . 85259E+07 -.16183E+06 -.89957E+05 O.OOOOOE+OO -.50603E+05 -.86194E+05 

4 0 . 85259E+07 0.16182E+06 0.89961E+05 0.00000E+00 0.49997E+05 O.OOOOOE+OO 
- . 85229E+07 -.16185E+06 -.90253E+05 -.18414E+04 -.50605E+05 -.86195E+05 



126 


The microstresses/strains are computed at the Gaussian center of all the elements. The first section 
of the file contains the microstresses at all the element-center points. The point location is identified 
by the X, Y and Z columns and the material property at the point is given by the ’yam no.’ column. 
A yam number of zero implies that the point is contained in the matrix material. The stresses/strains 
for a yarn point are listed in the yam LTS coordinate system (Fig. 4). However, the stresses/strains 
for a matrix point are listed in the global XYZ coordinate system. Therefore ’sigma LL’ for a yam 
point would mean the longitudinal stress component; and the normal stress component along the 
x-direction for a matrix point. The last section of the file lists the maximum tensile and compressive 
stresses/strains in the different constituent materials. 

4.1.3 Interpreting yarn. dat 

The element stiffness matrix is computed as the averaged stiffness of the constituent 
materials in that element (Appendix A. 1). For this purpose, it is required to determine the material 
property of the element Gauss points. The material property of a Gauss point is obtained in the 
following fashion. The yam volume is discretized into eight-node hexahedral isoparametric finite 
elements as shown in Fig. 7. The polygon points (used to reperesent the yarn cross-section in section 
2.1.3) are the nodes on the circumference of the yam. For a given cross-section, the center node 
coordinates are obtained by averaging the circumferential nodal coordinates. The finite elements are 
shaped as triangular prisms (a face of the hexahedron is collapsed to an edge). The length of the prism 
is governed by the spacing of the axial points, which represent the yarn path. An iterative algorithm 
is used to determine whether the point (whose material property is to be determined) is contained 
within the volume of a given element. If not, the procedure is repeated for all the elements in the 
yam. If the point does not belong to any of the yams, it is designated as a matrix point. If the point 
belongs to a yam element, the local yam direction is computed as the direction of the line along the 
length of the prismatic element. 



Figure 7. Finite element mesh for a yam 
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File yam.dat lists the finite element mesh information for all the yams. The data is given as 
the node number and the corresponding nodal coordinates; and then the element number and the 
connectivity for the element. The information is repeated for all the yams in serial order. The 
information is self-explanatory (since the data is printed out with appropriate headers), and hence 
is not explained in this manual. The user may transfer the data to a suitable graphics software 
package for visualizing the finite element mesh of the yams. 

Caution; The yam finite element mesh should not be confused with the finite element mesh of the 
unit-cell. The only purpose for generating the yam mesh is to determine the material property code 
at a given point (X,Y,Z). 

4.2 / iTEx-20 Output 

The output files ioxiiTE%— 20 are sam.out and yam.dat (see Appendix A.3). The format of 
sam.out is similar to that of textile. out’, and yam.dat is identical to that explained in section 4.1.3. 
Since the algorithm used to implement^ TE%—20 is based on an approximate analysis, the resulting 
stiffness matrix will not be symmetric. For the continuum problem (when ’solid’ option is used), 
the effective elastic constants are obtained from the stiffness matrix using the following two 
methods. In the first method the stiffness matrix is symmetrized, inverted to give the compliance 
matrix from which the elastic constants are predicted. In the second method, the (unsymmetric) 
stiffness matrix is inverted and the resulting compliance matrix is symmetrized. The results from 
both methods are printed in sam.out. 
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Appendix 

A.1 Computation of Inhomogenous Element Stiffness 

The element stiffness matrix (K e ) in /xTE%-10 is formulated as: 

K e = B t CB dV e 
J ye 

+1 +1 +1 

B t CB 171 d^drjdQ 

-1 -1 -1 
N N N 

= XX X WiWjWt B t C(Z, V ,0 J5I7I 

1=1 j=lk=l 

where V s is the domain of the element, B is the strain-displacement transformation matrix, C is the 
elasticity matrix, N is the number of Gauss points used for integration, W is the Gaussian integration 
weight factor and 171 is the determinant of the Jacobian. The material property at each Gauss point 
is determined, and the corresponding elasticity matrix is used to perform the volume 
integration. The element stiffness matrix thus represents the averaged properties of the constituent 
materials in that element. 




s 


A.2 File Flow Chart for fiTE%-10 


INPUT 


textile.inp 


yarn.inp 


luTEx-10 

(microtech 10) 


OUTPUT 


textile.out 


yarn.dat 


’stress’ or ’strain’ option I ’tstress’ or ’tstrain’ option 


stressl.dat stress6.dat 


stressl.dat stress6.dat , tstress.dat 


strainl.dat 7 \strainl.dat,..., strain6.dat, tstrain.dat 
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A.3 File Flow Chart for /uTEx-20 
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A.4 Stiffness Parameters in paramlO.var 

Required value for the global stiffness parameter NSTIFF in file paramlO.var. NX, NY and NZ 
denote the number of elements chosen to discretize the unit-cell in the three directions. 


NX 

NY 

NZ 

NSTIFF 

SOLID PLATE 

8 

8 

8 

549051 

373371 

8 

8 

9 

631587 

419043 

8 

8 

10 

714123 

464715 

8 

8 

11 

796659 

510387 

8 

8 

12 

879195 

556059 

8 

9 

8 

686883 

463512 

8 

9 

9 

790320 

520293 

8 

9 

10 

893757 

577074 

8 

9 

11 

997194 

633855 

8 

9 

12 

1100631 

690636 

8 

10 

8 

839691 

562869 

8 

10 

9 

966333 

631911 

8 

10 

10 

1092975 

700953 

8 

10 

11 

1219617 

769995 

8 

10 

12 

1346259 

839037 

8 

11 

8 

1007475 

671442 

8 

11 

9 

1159626 

753897 

8 

11 

10 

1311777 

836352 

8 

11 

11 

1463928 

918807 

8 

11 

12 

1616079 

1001262 

8 

12 

8 

1190235 

789231 

8 

12 

9 

1370199 

886251 

8 

12 

10 

1550163 

983271 

8 

12 

11 

1730127 

1080291 

8 

12 

12 

1910091 

1177311 

9 

8 

8 

693714 

471234 

9 

8 

9 

798015 

528879 

9 

8 

10 

902316 

586524 

9 

8 

11 

1006617 

644169 

9 

8 

12 

1110918 

701814 

9 

9 

8 

867999 

585147 

9 

9 

9 

998733 

656832 

9 

9 

10 

1129467 

728517 

9 

9 

11 

1260201 

800202 

9 

9 

12 

1390935 

871887 

9 

10 

8 

1061238 

710724 
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NX 

NY 

NZ 

NSTIFF 

SOLID PLATE 

9 

10 

9 

1221321 

797907 

9 

10 

10 

1381404 

885090 

9 

10 

11 

1541487 

972273 

9 

10 

12 

1701570 

1059456 

9 

11 

8 

1273431 

847965 

9 

11 

9 

1465779 

952104 

9 

11 

10 

1658127 

1056243 

9 

11 

11 

1850475 

1160382 

9 

11 

12 

2042823 

1264521 

9 

12 

8 

1504578 

996870 

9 

12 

9 

1732107 

1119423 

9 

12 

10 

1959636 

1241976 

9 

12 

11 

2187165 

1364529 

9 

12 

12 

2414694 

1487082 

10 

8 

8 

855225 

580419 

10 

8 

9 

983829 

651423 

10 

8 

10 

1112433 

722427 

10 

8 

11 

1241037 

793431 

10 

8 

12 

1369641 

864435 

10 

9 

8 

1070229 

720876 

10 

9 

9 

1231446 

809193 

10 

9 

10 

1392663 

897510 

10 

9 

11 

1553880 

985827 

10 

9 

12 

1715097 

1074144 

10 

10 

8 

1308633 

875733 

10 

10 

9 

1506063 

983163 

10 

10 

10 

1703493 

1090593 

10 

10 

11 

1900923 

1198023 

10 

10 

12 

2098353 

1305453 

10 

11 

8 

1570437 

1044990 

10 

11 

9 

1807680 

1173333 

10 

11 

10 

2044923 

1301676 

10 

11 

11 

2282166 

1430019 

10 

11 

12 

2519409 

1558362 

10 

12 

8 

1855641 

1228647 

10 

12 

9 

2136297 

1379703 

10 

12 

10 

2416953 

1530759 

10 

12 

11 

2697609 

1681815 

10 

12 

12 

2978265 

1832871 

11 

8 

8 

1033584 

700926 

11 

8 

9 

1189029 

786675 
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NX NY NZ NSTEFF 

SOLID PLATE 


11 

8 

10 

1344474 

872424 

11 

8 

11 

1499919 

958173 

11 

8 

12 

1655364 

1043922 

11 

9 

8 

1293573 

870699 

11 

9 

9 

1488459 

977376 

11 

9 

10 

1683345 

1084053 

11 

9 

11 

1878231 

1190730 

11 

9 

12 

2073117 

1297407 

11 

10 

8 

1581876 

1057896 

11 

10 

9 

1820559 

1187679 

11 

10 

10 

2059242 

1317462 

11 

10 

11 

2297925 

1447245 

11 

10 

12 

2536608 

1577028 

11 

11 

8 

1898493 

1262517 

11 

11 

9 

2185329 

1417584 

11 

11 

10 

2472165 

1572651 

11 

11 

11 

2759001 

1727718 

11 

11 

12 

3045837 

1882785 

11 

12 

8 

2243424 

1484562 

11 

12 

9 

2582769 

1667091 

11 

12 
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